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ABSTRACT 

System  identification  concerns  the  mathematical, modeling  of  a  system  based  upon 
its  input  and  output.  It  allows  the  development  of  a  mathematical  description  when  all 
that  is  available  is  the  result  of  a  process  or  the  output  of  a  system  and  not  the  process 
or  system  itself. 

The  purpose  of  this  thesis  is  to  develop  algorithms  for  modeling  systems  as 
autoregressive-moving-average  processes  using  the  method  of  instrumental  variables,  a 
modification  of  the  ordinary  least-squares  technique,  and  a  multichannel  method  based 
upon  processing  the  input  and  output  data  by  separate  infinite-impulse-response  filters. 
The  methods  developed  are  tested  by  computer  simulation  using  several  second  and 
third-order  test  cases  and  the  results  are  presented. 
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I.     INTRODUCTION 

A.     SYSTEM  IDENTIFICATION  BASICS 

System  identification  concerns  the  modeling  of  systems  as  sets  of  mathematical 
equations  based  upon  the  input  and  output  of  the  system.  [Ref.  1:  pp.  3-6].  It  allows 
a  model  to  be  developed  when  all  that  is  available  is  the  result  of  a  process  or  the  output 
of  a  system  and  not  the  process  or  the  system  itself.  System  identification  is  an  impor- 
tant area  of  study.  Solution  of  the  modeling  problem  offers  many  alternatives  for  the 
continued  study  of  the  system.   Among  these  are: 

•  Nondestructive  analysis  of  the  system. 

•  Simulation  studies  using  the  model. 

•  Easy  adaptation  of  the  model  to  changing  system  environment. 

•  Spectral  analysis  of  the  system. 

Modeling  can  simulate  the  system's  operation  at  a  fraction  of  the  cost  of  actual 
system  operation.  Complex  operations  not  possible  with  the  actual  system  for  fear  of 
damaging  it  or  personal  injury  can  be  simulated.  This  can  expose  how  the  system  will 
operate  in  adverse  conditions  not  normally  experienced.  In  speech  processing,  modeling 
the  speech  process  has  the  potential  for  significantly  reducing  the  amount  of  information 
necessary  to  store  in  order  to  reproduce  the  speech. 

The  modeling  process  shown  in  Figure  1  on  page  2  assumes  the  unknown  system's 
input  and  the  output  data  are  available  for  processing.  In  many  cases,  if  the  system's 
input  is  unknown  or  data  is  not  available,  a  white  noise  input  can  be  used  in  its  place. 
The  modeling  process  uses  the  input  and  output  data  to  find  a  set  of  parameters  which 
closely  approximate  the  operation  of  the  system.  The  better  the  identification  technique, 
the  more  closely  the  model  follows  the  performance  of  the  actual  system. 

Many  types  of  models  are  available.  This  thesis  investigates  a  linear  parametric 
model  that  can  be  described  by  difference  equations.  This  type  of  model  lends  itself  well 
to  simulation  on  a  digital  computer.  The  frequency  characteristics  of  the  system  deter- 
mined from  the  parameters  of  these  types  of  models  are  more  accurate  than  what  can 
be  determined  from  classical  means  such  as  FFTs.  This  is  because  classical  methods  use 
windows  which  assume  data  beyond  their  extent  is  zero  [Ref.  2:  p.  173].  This  is  not  a 
realistic  assumption.    Models  in  this  category  include  the  moving-average  (MA)  model, 
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Figure  1.      System  identification  problem 

the  autoregressive  (AR)  model,  and  the  autoregressive-moving-average  (ARMA)  model. 
In  the  frequency  domain,  MA  processes  are  characterized  by  sharp  nulls  and  smooth 
peaks  and  AR  processes  are  characterized  by  smooth  nulls  and  sharp  peaks.  ARMA 
processes  have  sharp  peaks  and  sharp  nulls  [Ref.  2:  p.  173].  An  advantage  of  the  MA 
process  is  its  inherent  stability.  An  advantage  of  the  AR  process  is  the  large  number  of 
algorithms  already  available  for  modeling  systems.  An  advantage  of  the  ARMA  process 
is  that  it  uses  far  fewer  parameters  than  either  the  MA  or  AR  process  alone  to  model  a 
system.   This  satisfies  the  general  requirement  to  reduce  the  complexity  of  the  model. 

In  addition  to  a  large  variety  of  models,  there  are  two  processing  modes:  block  and 
sequential. 

Block  processing  uses  a  fixed  length  block  of  data  in  the  parameter  estimation 
process.  It  ignores  data  before  and  after  the  block.  This  is  not  a  real  time  processing 
method  because  all  data  must  be  available  before  processing  can  start.  Block  processing 
generally  involves  inversions  of  data  matrices  whose  sizes  are  on  the  order  of 
(iV+  M)  x(.Y+  M)  where  N  is  the  order  of  the  AR  process  and  M  is  the  order  of  the 
MA  process. 


Sequential  processing  uses  new  data  to  update  the  parameter  estimations.  It  starts 
by  initializing  an  estimate  of  the  inverse  of  the  data  covariance  as  as  a  diagonal  matrix. 
It  uses  each  new  data  point  to  update  this  matrix.  Then  it  updates  the  parameter  esti- 
mates using  the  updated  inverse  data  covariance  matrix.  It  is  a  real  time  method.  The 
algorithm  to  implement  the  sequential  processing  method  is  generally  more  complex 
than  the  block  method  but  less  computationally  intensive  because  the  matrix  inversions 
are  not  required. 

This  thesis  concerns  only  systems  represented  by  discrete  time  data  uniformly  sam- 
pled at  a  sufficient  rate  to  meet  the  Nyquist  criteria. 

The  work  in  this  thesis  assumes  that  the  input  data  is  a  wide-sense  stationary  ran- 
dom sequence.  Tests  of  the  algorithms  used  a  pseudorandom  Gaussian  input  with  a 
mean  of  zero  and  a  variance  of  one. 

B.  PROBLEM  STATEMENT 

The  purpose  of  this  thesis  is  to  develop  algorithms  for  modeling  systems  as  ARM  A 
processes  using  the  method  of  instrumental  variables  (IV)  and  a  multichannel  approach. 
Tests  of  the  methods  will  be  conducted  to  determine  the  accuracy  of  their  results  and  the 
speed  with  which  they  converge. 

The  IV  approach  is  a  modification  of  the  method  of  ordinary  least  squares.  This 
approach  is  developed  first  as  a  block  processing  case  and  then  converted  to  a  sequential 
processing  case.   Tests  are  conducted  of  only  the  sequential  processing  case. 

Using  a  multichannel  scheme  allows  the  input  and  output  data  of  the  unknown 
system  to  be  processed  separately.  This  reduces  the  sizes  of  the  data  matrices  involved 
in  the  modeling  process.  Both  block  and  sequential  processing  cases  are  formulated  but 
only  the  block  processing  case  is  tested. 

C.  OVERVIEW  OF  THESIS 

Chapter  2  is  about  ARMA  modeling.  It  also  presents  a  detailed  derivation  of  the 
method  of  ordinary  least  squares  because  it  forms  the  basis  on  which  other  modeling 
techniques  depend. 

Chapter  3  presents  a  modified  least-squares  approach  called  the  method  of  instru- 
mental variables.  It  is  attractive  due  to  its  simplicity  and  good  noise  performance. 
Chapter  3  presents  results  of  using  this  method  on  several  second  and  third-order  test 
svstems. 


Chapter  4  presents  a  new  multichannel  approach  to  ARM  A  modeling.  This  ap- 
proach is  presented  in  block  and  sequential  processing  forms.  This  chapter  also  presents 
several  adaptations  of  the  block  form  which  improve  its  speed  of  convergence. 

Chapter  5  contains  a  summary  of  the  thesis  and  lists  topics  for  further  research. 

The  appendix  contains  the  programs  used  to  test  the  sequential  IV  algorithm  and 
the  block  multichannel  iterative  algorithm.  Subroutines  common  to  both  programs  are 
grouped  together  and  listed  at  the  end  of  the  appendix. 


II.     ARMA  MODELING 

A.     ARMA  PROCESSES 

Modeling  as  an  autoregressive-moving-average  (ARMA)  process  has  the  potential 
for  achieving  a  close  fit  to  the  system  using  a  reduced  order  over  that  which  a  moving 
average  or  an  autoregressive  model  alone  could  achieve.  ARMA  modeling  is  concerned 
with  finding  a  set  of  AR  parameters  and  MA  parameters  which  combined  describe  an 
ARMA  process  that  approximates  the  characteristics  of  a  target  system. 

The  general  form  of  the  ARMA  model  is  shown  in  Figure  2  on  page  6.  The  output 
at  time  n.  y(n),  is  a  linear  combination  of  past  outputs  and  past  and  present  inputs. 
The  a,  and  b:  are  constants  referred  to  as  tap  weights.  The  at  parameters  form  the  MA 
part  of  the  ARMA  model.  The  b,  parameters  form  the  AR  part.  In  equation  form  the 
output  of  the  ARMA  system  is  represented  by  the  following  difference  equation: 

iV  M 

(=  1  (=0 

where  N  is  the  order  of  the  AR  part  of  the  ARMA  model  and  M  is  the  order  of  the  MA 
part  of  the  ARMA  model.  This  means  the  ARMA  output  at  the  current  time  depends 
on  the  last  ;V  values  of  the  ARMA  output.  The  A'  b,  weighting  parameters  determine 
exactly  how  the  new  output  depends  on  the  past  outputs.  The  M  a,  weighting  parame- 
ters determine  how  the  new  output  depends  on  the  current  and  M  —  1  past  inputs. 
Equation  (2.1)  in  vector  form  becomes: 

y  =  xT6  (2.2) 

where  x  is  a  (Ar+  M  +  1)  x  1  vector  of  input  and  output  data  values  given  by: 

x  =  l-y(n-\)     -y{n-2)     ...      -y(n  -  N)   x(n)   x(n  -  1)     ...     *(«-  M)f  (2.3) 

and  6  is  a  (X  +  M  +  1)  x  1  vector  of  the  AR  and  MA  tap  weights  given  by: 

e  =  [bx    b2     ...     bN   a0    a,     ...     a„f  (2.4) 
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Figure  2.      Structure  of  the  ARMA  model 

For  N  +  L  —  I  data  points  available  for  y  and  M  +  L  data  points  available  for  u,  we 
can  write  a  block  equation  which  gives  the  value  of  the  output  at  progressive  sampling 
times: 


y(n-L+\) 
y(n-L  +  2) 


A") 


\T(n+  1), 
xT(n  +  2) 


x\n  +  L) 


bN 


lM 


(2.5) 


The  j01  row  in  equation  (2.5)  is  the  value  ofj  at  time  n  +  i  based  on  output  data  available 
through  time  n  +  i  —  \  and  input  data  available  through  n  +  i.  The  Ith  row  is  identically 
equation  (2.2)  at  time  n  +  i.    In  vector  form  equation  (2.5)  becomes: 

y  =  X0  (2.6) 

where  6  is  defined  in  equation  (2.4);  y,  the  vector  of  output  values,  is  given  by: 

y  =  [y(«-L+l)  A"-L  +  2)     ...    y(n)f  (2.7) 


and  X  is  a  partitioned  matrix  with  rows  comprised  of  data  vectors  exactly  like  equation 
(2.3)  only  shifted  in  time.  At  successive  sampling  times,  when  new  data  is  obtained,  data 
used  to  calculate  the  previous  output  shifts  one  column  to  the  right.  The  new  data  fills 
in  the  left  most  y  and  u  columns.   The  matrix  X  is  given  by: 


X  = 


-y{n  -  L) 
-y(n-L+\) 


-y(n-l) 


-y{n-rj+\)      u{n  -  L  +  1) 
-y(n  -  rj  +2)      u(n  -  L  +  2) 


-M«-iV) 
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(2.8) 


where  r\  is  defined  as  AT+  L  and  /x  is  defined  as  M  +  L. 


If  the  a,  and  b,  are  estimates  of  the  true  values  of  the  AR  and  MA  parameters,  then 
the  filter  output  will  be  an  estimate  of  the  true  output.  We  use  a  hat  over  a  variable  (for 
example,  y)  to  indicate  an  estimated  value.  Rewriting  equation  (2.6)  using  the  estimated 
ARM  A  parameters  results  in: 


y  =  X0  (2.9) 


where  0  is  defined  as: 


A  A  A 


0  =  l>i    b2     ...     bN   aQ    a,     ...     aM]r  (2.10) 

and  y  is  now  the  vector  of  estimated  output  values  and  is  given  by: 

y=[y(n-L+l)  y{n-L  +  2)     ...    y(n)f  (2.11) 

Up  until  now,  we  have  discussed  estimating  the  output  of  a  system  given  its  input, 
past  output,  and  an  estimate  of  the  parameters  which  describe  it.  If,  however,  we  know 
the  output  and  input  of  the  system,  based  on  these  equations,  we  can  use  them  to  gen- 
erate a  set  of  a,  and  b,  which  produces  an  ARMA  output  which  is  the  best  possible  es- 
timate of  the  system  output.  Then  the  a,  and  b,  will  be  optimal  parameters  for  describing 
the  operation  of  the  unknown  system  as  an  ARMA  process. 

B.     METHOD  OF  ORDINARY  LEAST-SQUARES 

In  this  thesis  we  use  the  method  of  ordinary  least-  squares  as  the  means  of  finding 
the  optimum  set  of  ARMA  parameters.  It  is  a  well  known  modeling  technique.  It  offers 
the  advantage  of  being  widely  used  in  the  scientific  community  for  a  variety  of  modeling 
problems.  It  has  been  applied  successfully  to  a  large  number  of  modeling  problems  with 
good  results  and  has  been  successfully  applied  to  classes  of  problems  for  which  other 
methods  have  failed.   [Ref.  3:  p.  4] 

To  apply  the  method  of  ordinary  least-squares  to  system  identification  we  form  the 
error  between  the  actual  system  output  and  the  estimated  output  generated  by  the 
ARMA  model.   This  error  is  given  by: 

s  =  y-y  =  y-X0  (2.12) 

where  y  is  the  vector  of  the  actual  system  outputs  given  by  equation  (2.7)  and  y  is  the 
vector  of  ARMA  outputs  given  by  equation  (2.11).   [Ref.  1:  p.  176] 


We  then  let  the  sum  of  the  squares  of  the  errors  at  the  instances  of  time  the  meas- 
urements of  the  data  were  taken  become  a  measure  of  how  well  the  estimates  approxi- 
mate the  true  system  outputs.  This  measure  of  performance,  or  cost  function,  is  denoted 
J  .    It  is  written  in  equation  form  as: 


n+L 

J  =    >    e]  =  &t&  (2.13) 

i=n+\ 

Replacing  the  error  in  equation  (2.13)  with  its  equivalent  expression  from  equation 
(2.12)  results  in: 

J  =  yry  +  6TXXr6  -  20rXy  (2.14) 

Equation  (2.14)  shows  that  the  performance  measure  is  a  function  of  the  estimated  pa- 
rameters. The  criterion  is  to  minimize  the  measure  of  performance  by  taking  its  deriva- 
tive with  respect  to  the  parameter  estimates  and  setting  it  equal  to  zero.  Then  equation 
(2.14)  becomes: 

SL  =  o  =  o  +  2XXT6  -  2Xry  (2.15) 

60 


Solving  for  6,  the  parameters,  gives  us  the  result: 

d  =  (XTX)~]XTy  (2.16) 


Equation  (2.16)  is  the  ordinary  least-squares  solution  for  the  optimum  ARM  A  pa- 
rameters. It  provides  the  best  possible  description,  in  a  least-squares  sense,  of  the  data 
source.  The  resulting  parameters  provide  the  closest  fit  to  the  actual  input  and  output 
data  of  the  system  in  the  sense  of  least-squares  errors. 

Equation  (2.16)  uses  a  block  processing  approach.  The  product  of  X'X  must  be 
formed  and  then  inverted  in  order  to  calculate  0  .  In  addition  to  being  computationally 
intensive,  the  estimate  cannot  be  updated  when  new  data  becomes  available  without  re- 
calculating (X^)-1  .  A  sequential  update  which  does  not  require  (XOC)-1  to  be  recalcu- 
lated is  presented  in  the  next  chapter  in  the  context  of  the  instrumental  variable  method 
of  least-squares. 


III.     INSTRUMENTAL  VARIABLE  METHOD  OF  SYSTEM 

IDENTIFICATION 

A.     INTRODUCTION 

The  instrumental  variable  (IV)  method  of  system  identification  is  a  variation  of  the 
method  of  ordinary  least-squares.  Its  attraction  over  ordinary  least-squares  is  that  there 
is  no  bias  in  estimating  the  parameters  when  dealing  with  noise  [Ref.  4:  p.  406].  Also, 
this  method  is  known  to  yield  consistent  estimates  and  remains  as  easy  to  use  as  the 
method  of  ordinary  least-squares  [Ref.  3:  p.  119]. 

When  an  additive  noise  term  is  present  in  the  observable  output,  y(n),  the  output  is 
given  by: 


e 


y(n)  =  w(n)  +  V(n)  (3.1) 

Here  w(n)  represents  the  actual  output  of  the  system  and  v(n)  represents  the  noise.  When 
this  noise  has  a  non-zero  mean,  using  the  noise  corrupted  output  to  model  the  unknown 
system  by  the  ordinary  least-squares  approach  leads  to  inaccurate  estimates  of  its  pa- 
rameters. The  parameters  are  referred  to  as  biased  estimates.  [Ref.  3:  p.  119,  Ref.  1:  pp. 
192-193,  and  Ref.  5:  p.  704]. 

The  IV  method  shown  in  Figure  3  on  page  11  generates  an  estimate  of  the  un- 
known system's  output  by  processing  the  input  data  through  an  auxiliary  model  which 
closely  approximates  the  unknown  system.  In  our  implementation  of  the  IV  method, 
the  auxiliary  model  is  an  ARMA  model.  Its  output  is  free  of  the  noise  affecting  the 
unknown  system.  The  IV  method  uses  the  auxiliary  model  output  (estimate),  w,  to  cal- 
culate the  parameters  of  the  unknown  system.  Therefore  the  IV  parameter  estimates  are 
not  biased  like  those  generated  by  the  method  of  ordinary  least-squares. 

The  IV  method  assumes  the  existence  of  a  matrix  Z  composed  of  the  auxiliary 
model's  input  and  output  data  which  has  the  following  two  properties  [Ref.  4:  p.  406]: 

lim  -—  TJz  =  0  (3.2) 

A'-»oo    A 


lim-frZrX  =  Q  (3.3) 

where  e  is  the  error  in  fitting  the  parameter  estimates  to  the  data  and  is  given  by: 
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v(n) 


u(n) 


Ls 


SYSTEM 
A(z)/B(z) 

w(n)     IT      y( 

Ci) 

AUXILIARY  MODEL 

A            A 

A(z)/B(z) 

w(n) 

1 

Figure  3.      Modeling  by  the  instrumental  variable  method 


E  =  V  -  X6 


/K 


(3.4) 


and  Q  is  a  nonsingular  square  matrix. 

The  first  property  means  Z  is  orthogonal  to  the  error.  This  leads  to  the  cancellation 
of  the  bias  term  inherent  in  ordinary  least-squares  techniques  [Ref.  4:  p.  406].  The  sec- 
ond property  ensures  the  inverse  of  Z^  exists.  Z  is  assumed  to  have  the  same  structure 
and  size  as  the  data  matrix  X  in  equation  (2.8).  Its  contents  differ  in  that  the  noise 
corrupted  system  output  y(n)  in  X  is  replaced  by  the  output  of  the  auxiliary  model  \v(n) 
in  Z.   The  new  data  matrix  Z  is  given  by: 


Z  = 


—  w(n  —  L) 
-  w(n  -L+l) 


-w{n-  1) 


-  w{n-r}  +1)      u{n-  L+  1) 

-  w{n  -  rj  +2)      u{n-  L  +  2) 


—  w(n  —  iY) 


u(n) 


u(n  - 

-n+1) 

u(n  - 

■M+2) 

u(n 

-\f) 

(3.5) 
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where  r\  is  defined  as  .V  4-  L  and  fx  is  defined  as  M  +  L.  Comparing  X  in  equation  (2.8) 
and  Z  in  equation  (3.5),  we  note  the  substitution  of  w(n)  for  y{n).  Thus,  we  are  now 
using  estimates  of  the  true  output  w{n)  instead  of  noise  corrupted  samples  y(n). 

To  incorporate  Z  into  the  parameter  estimation  process  we  begin  with  equation 
(2.12),  which  we  rewrite  as: 

y  =  X0  +  s  (3.6) 

This  equation  says  that  the  estimates  of  the  output,  given  by  X0,  differ  from  the  actual 
outputs,  y,  by  some  fitting  error  s.    Multiplying  equation  (3.6)  by  ZT  yields: 

Zry  =  ZrX0  +  Zr£  (3.7) 

Equation  (3.3)  ensures  that  Z'X  can  be  inverted.    Solving  for  6  results  in: 

B  =  (ZTX)~lZTy-{ZTX)~lZTs  (3.8) 

The  (Z'X^Z7}7  term  in  equation  (3.8)  is  the  IV  estimate  of  the  parameters.  It  is  written 
as: 

0IV=(ZTX)-]ZTy  (3.9) 

The  {Z1XyxZT&  term  in  equation  (3.8)  represents  a  potential  bias  in  the  estimate.  The 
first  property  of  the  Z  matrix,  given  in  equation  (3.2),  ensures  this  bias  goes  to  zero, 
asymptotically.   Applying  this  property,  equation  (3.8)  can  be  rewritten  as: 

8  =  (Z7X)-1Z7y-8/F  (3.10) 

Equation  (3.10)  gives  an  unbiased  estimate  of  the  ARMA  parameters.  [Ref.  1:  pp. 
192-193]: 

Other  least-squares  methods  avoid  the  bias  inherent  in  ordinary  least-squares  but 
they  are  more  complicated  than  the  IV  method  to  implement  [Ref.  3:  p.  119].  Although 
this  thesis  does  not  attempt  an  analysis  of  the  IV  method  in  the  presence  of  noise,  any 
practical  system  identification  technique  must  deal  with  noise.  Hence,  the  attraction  of 
and  the  desire  to  use  the  IV  method. 

Equation  (3.10)  represents  the  block  processing  case.  It  assumes  /V+  L  —  1  output 
samples  and  M  +  L  input  samples  are  available.   These  samples  are  used  to  calculate  an 
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estimate  of  the  parameters.  Samples  beyond  this  range  are  not  included  in  the  esti- 
mation process.  Block  processing  involves  multiplication  of  two  L  x  (A"  +  M  +  1)  ma- 
trices to  form  a  third  matrix.  Then  this  third  matrix  must  be  inverted.  This  is  a 
computationally  intensive  process.  In  what  follows,  we  present  a  sequential  algorithm 
to  compute  6JV  which  avoids  matrix  inversions. 

B.     SEQUENTIAL  LEAST-SQUARES  ESTIMATION  USING  INSTRUMENTAL 
VARIABLES 

A  sequential  process  for  estimating  the  parameters  of  an  unknown  system  requires 
fewer  computations  than  a  block  process.  In  a  manner  similar  to  that  presented  in  Hsia 
[Ref.  3:  pp.  22-25]  for  the  general  least-squares  case,  the  block  IV  estimation  process 
described  above  can  be  converted  into  a  sequential  IV  estimation  process.  Using  the 
sequential  process  also  allows  the  coefficients  to  be  updated  based  on  the  new  data  that 
becomes  available. 

The  derivation  of  the  sequential  estimation  procedure  consists  of  two  parts.  The 
first  part  is  the  derivation  of  an  equation  to  update  the  data  matrix,  Q(m  +  1) ,  based 
on  the  previous  data  matrix,  Q(m) ,  and  the  new  data:  w{m) ,  y(m) ,  and  u(m  +  1) 
where  m  represents  the  iteration.  The  second  part  involves  developing  an  equation  for 
updating  the  estimate  of  the  parameters,  0/K(m  +  1),  based  on  the  previous  estimate, 
6n(m) ,  the  previous  data  matrix,  Q(m) ,  and  the  new  data:  w{m) ,  y(m) ,  and 
u(m  +  1)  . 

Define  the  data  matrix  Q(m)  to  be: 

Q(m)  =  [ZXr1  (3.11) 

where  Zm  is  given  by  equation  (3.5)  and  Xm  is  given  by  equation  (2.8).  The  property  of 
equation  (3.3)  assures  that  Q  exists.  Note  that  Q(m)  includes  output  data  available 
through  m  and  input  data  available  through  m  +  1  .  Since  both  Zm  and  Xm  are 
m  x  (V  +  \f)  matrices,  Q(m)  will  be  a  (Ar  +  M)  x  (N  +  Af)  matrix.  As  the  number  of  rows 
of  Z  and  X  increase  to  accommodate  the  increasing  numbers  of  data  points,  the  size  of 
Q  will  remain  the  same.  At  the  next  sample  time,  i.e.,  at  m  +  1,  the  data  matrix  becomes: 

Q(m  +  1)  =  [z£+1Xm+1]-]  (3.12) 

where  the  data  matrices  at  m  +  1  are  given  by: 


13 


7 

*-'m 

7         = 

**m+l 

zT{m+  1) 

^m 

*-m+l  = 

.... 

xT{m+  1) 

(3.13) 


(3.14) 


and  zT(m  +  1)  and  xT(m  +  1)  are  vectors  which  contain  the  most  recent  data  values. 
Substituting  equations  (3.13)  and  (3.14)  into  equation  (3.12)  and  expanding,  results  in: 


Q(m+  1)  = 


[Z;     z(m+l)] 


X, 


xy(m+  1) 


-i 


(3.15) 


Expanding  further  yields: 


Q(m  +  1)  =  \ZTmXm  +  z(m  +  l)xr(m  +  1)] 


-l 


(3.16) 


In  equation  (3.16)  we  see  that  two  terms  make  up  the  new  data  matrix.  The  Z£Xm  term 
is  all  the  data  that  was  available  through  time  m.  The  z(m  +  \)xT(m  +  1)  term  contains 
the  new  data.  To  perform  the  inversion,  let  A  =  Z£Xm,  B  =  z(m+1),  C=l  and 
D  =  xT(m  +  1)  .   Then  by  the  matrix  inversion  lemma: 


Q(m  +  1)  =  A"1  -  A_1B(C_1  +  DA"lB)~1DA"1 


(3.17) 


Substituting  the  appropriate  expressions  for  A,  B,  C,  and  D  back  into  equation  (3.17) 
yields  the  equation: 


Q(m  +  1)  -  (Z£Xm)-]  -  (ZX)"'z(m  +  1) 

.    [l+xT(m+l)(ZlXm)-'z(m+l)Y] 

.  xV+ixzXr1 


(3.18) 


Substituting  Q(m)  for  (Z^XJ-1  reduces  equation  (3.18)  to: 

Q(m  +  1)  =  Q(m)  -  Q(m)z(m  +  1)[1  +  xT{m  +  l)Q(w)z(w  4-  1)]_1 
.    xT{m+  l)Q(m) 


(3.19) 
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This  completes  the  first  step  of  the  derivation.  Equation  (3.19)  expresses  Q  at  time 
m  +  1  in  terms  of  the  old  Q  and  the  new  data.  The  term  in  the  brackets  is  a  scalar. 
Computational  intensity  has  been  reduced  because  a  large  matrix  does  not  have  to  be 
generated  and  its  inverse  does  not  have  to  be  calculated. 

Continuing  with  the  derivation,  the  estimate  BIV  for  data  available  through  m  can 
be  written  as: 


0lv(m)  =  (Z„?X  J    Zmym 

A 

The  estimate  0IV  for  data  available  through  m  +  1  can  be  written  as: 


(3.20) 


Bn{m  +  1)  —  {Zm+]Xm+])    Zw+1yWJ+1 


(3.21) 


Substituting  equation  (3.12)  into  equation  (3.21)  results  in  an  expression  for  the  estimate 
of  the  parameters  in  terms  of  the  new  data  matrix  and  all  the  available  data  given  by: 


0lv(m+\)  =  Q(m+\)Zll+tfm+l 


(3.22) 


yn+l)  =  Q(ffl+l)[Z;     z(m+l)] 


_y(m+  1)_ 


(3.23) 


eiv{m  +  1)  =  Q(m  +  l)[Z;ym  +  z(m  +  l)y(m  +  1)] 


(3-24) 


Substituting  for  Q(m  +  1)  from  equation  (3.19)  and  expanding  results  in: 


0ltAm  +  1)  =  Q(m)Z!nym 

-  Q(m)z(m  +  1)[1  +  xT(m  +  l)Q(m)z(m  +  l)TlxT(m  +  l)Q(m)Z^ym 


+  Q(m)z(m+  \)y{m+  1) 

-  Q(m)z(m  +  1)[1  +  xT(m  +  l)Q(m)z(m  +  l)]"1 

.    xT(m  +  l)Q(m)z(m  +  \)y{m  +  1) 


(3.25) 


Although  somewhat  lengthy,  this  equation  has  the  desired  form.    To  simplify  it,  its  last 
two  terms  can  be  arranged  into  the  form: 


i-l_.7V 


Q(m)z(m  +  1){1  -  [1  +  x'(m  +  l)Q(m)z(m  +  1)]    \'{m  +  l)Q(m)z(m  +  1)}      3  2„ 

•  y{m+  1) 
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The  terms  within  the  braces  can  be  thought  of  as  the  result  of  a  previous  application  of 
the  matrix  inversion  lemma  with  A~'=l,  B=l,  C_1  =  1  and 
D  =  xr(m  +  l)Q(m)z(m  +  1)  .    Reversing  the  lemma  results  in: 

Q(m)z{m  +  1)[1  4-  x\m  +  [)Q(m)z{m  +  1  )]">■(>"  +  1)  (3.27) 

Replacing  the  last  two  terms  in  equation  (3.25)  with  this  result  gives  us: 

0n(m  +  1)  =  Q{m)Zlym  -  Q(m)z(m  +  1)[1  +  xT(m  +  l)Q(m)z(m  +  1)]"' 

.    x\m  +  l)Q(m)Z£ym  +  Q(m)z(m  +  1)  (3.28) 

.    [1  +  x\m  +  l)Q(m)z(m  +  1)]"V('»  +  1) 

Factoring  Q(m)z{m  +  1)  and  [1  +  xT(m  +  l)Q(m)z(m  +  1)]_1  from  the  last  two  terms  re- 
duces equation  (3.28)  further  to: 

6n(m  +  1)  =  Q(«)Zjyjn  +  Q(m)z(m  +  1)[1  +  xT(m  +  \)Q(m)z(m  +  l)]"1 
.    0(m  +  1)  -  xr(m  +  l)Q(m)Zjym] 

Substituting  equation  (3.11)  into  equation  (3.20)  and  then  equation  (3.20)  into  equation 
(3.29)  yields  the  final  form  for  the  update  of  the  estimate  of  the  parameters: 

Blv{m  +  1)  =  0n{m)  +  Q(,n)z(m  +  1) 

a  (3.30) 

.    [I  +  xT{m  +  l)Q(m)z(m  +  1)]" x\y{m  +  1)  -  x\m  +  1)9/K(m)] 

This  is  the  desired  result  for  updating  the  estimate  of  the  parameters.  Note  that  like 
equation  (3.19),  the  matrix  inversion  of  equation  (3.21)  has  been  reduced  to  inversion 
of  a  scalar.  Equation  (3.30)  describes  the  update  of  dIV(m+  1)  in  terms  of  the  previous 
estimate  of  the  parameters,  6n{/n).  the  previous  data  matrix,  Q(w),  and  the  new  data: 
w(m) ,    y(m) ,    and  u(m  +  1)  . 

C.     TESTING  THE  SEQUENTIAL  INSTRUMENTAL  VARIABLE  ALGORITHM 

Equations  (3.19)  and  (3.30)  above  comprise  the  sequential  IV  algorithm.  Several 
tests  of  this  algorithm  were  made  using  second  and  third-order  filters  as  unknown  sys- 
tems. Tests  were  run  via  computer  simulation  using  the  filters  to  generate  the  output 
data.  A  Gaussian  random  process  with  zero  mean  and  unit  variance  was  used  as  the 
input.    The  input  was  produced  by  IMSL  subroutine  GGXML.    Graphs  were  created 
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using  DISSPLA.    Table   1  on  page  17  shows  pole  and  zero  locations  as  well  as  numera- 
tor and  denominator  parameters  for  the  test  filters. 


Table   1.     TEST  SYSTEMS  FOR  THE  IV  MODELING  METHOD 


TEST 
FILTER 

LOCATIONS 
OF  POLES 

LOCATIONS 
OF  ZEROS 

AR 
PARAMETERS 

MA 
PARAMETERS 

T2 

0.445  + jO.228 

0.445  -  jO.228 

0.4+J1.273 

0.4  -  jl. 273 

1.0 
-0.89 
0.25 

0.5 
-0.4 
0.S9 

T2N 

0.445  +  J0.228 
0.445  -jO.228 

0.4  +  jO.S 
0.4  -jO.8 

1.0 
-0.89 
0.25 

1.0 
-0.80 
0.80 

T3 

0.6605 

0.6647 +  J0.502 
0.6647-J0.502 

-1.0 
-1.0 
-1.0 

1.0 
-1.99 

1.57 
-0.45S 

0.0154 
0.0462 
0.0462 
0.01 5  J 

Results  of  the  tests  are  shown  in  graphical  form  in  Figure  4  on  page  18,  Figure  5 
on  page  19.  and  Figure  6  on  page  20.  Dashed  lines  indicate  the  true  values  of  the  pa- 
rameters.   Solid  lines  are  the  IV  method's  estimates. 

For  both  second-order  test  cases  shown,  the  algorithm  converged  quickly  and 
produced  accurate  results.  For  the  third-order  test  case,  convergence  took  longer  but 
the  values  were  accurate.  A  third-order  system  is  more  complex  than  a  second-order 
system,  so  conceivably  it  would  require  more  iterations  to  converge.  The  number  of  it- 
erations required  is  of  the  same  order  as  the  method  of  ordinary  least-squares. 

Table  2  on  page  21  contains  the  IV  algorithm's  best  estimates  of  the  parameters  and 
the  number  of  iterations  required  to  converge  to  those  estimates.  It  also  shows  the  ab- 
solute and  percent  errors  from  the  true  parameters. 
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Figure  4.      Second-order  test  case  T2.  (A)  MA  parameters.   (B)  AR  parameters. 
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Figure  5. 


Second-order  test  case  T2N.  (A)  MA  parameters.  (B)  AR  parameters. 
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Figure  6.      Third-order  test  case  T3.  (A)  MA  parameters.   (B)  AR  parameters. 
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Table  2.     COEFFICIENT  ESTIMATES  BY  THE  IV  MODELING  METHOD. 


TEST 

PARAMETER 

ABSOLUTE 

PERCENT 

ITERATIONS 

FILTER 

ESTIMATE 

ERROR 

ERROR 

T2 

0.500 

0.0 

0.0 

10 

-0.396 

+  0.004 

0.10 

0.SS8 

-0.002 

'      0.22 

1.000 

0.0 

0.0 

-0.8S8 

+  0.002 

0.22 

0.244 

-0.006 

2.40 

T2N 

1.000 

0.0 

0.0 

10 

-0.794 

+  0.006 

0.750 

0.794 

-0.006 

0.750 

1.000 

0.0 

0.0 

-0.SS7 

+  0.003 

0.34 

0.243 

-0.007 

2.80 

T3 

0.0154 

0.0 

0.0 

1000 

0.0466 

+  0.0004 

0.87 

0.0475 

+  0.0013 

2.S1 

0.0169 

+  0.0015 

9.74 

1.000 

0.0 

0.0 

-1.96 

-0.03 

3.0 

1.532 

-0.040 

2.01 

0.4379 

-0.0204 

4.45 

IV.     SYSTEM  IDENTIFICATION  USING  AN  ITERATIVE 
MULTICHANNEL  APPROACH 

A.  INTRODUCTION 

This  chapter  presents  an  alternate  system  identification  method,  the  iterative  multi- 
channel approach.  This  approach  differs  from  the  IV  method  and  the  method  of  ordi- 
nary least-squares  presented  in  the  preceding  chapters  in  that  it  processes  the  input  and 
output  data  from  the  unknown  system  in  separate  channels.  In  its  block  processing 
form  one  advantage  over  the  IV  and  ordinary  least-squares  methods  is  a  reduction  in  the 
sizes  of  the  data  matrices.  As  a  result,  the  computational  complexity  of  the  multichannel 
algorithm  is  on  the  order  of  M1  +  A'2  ,  where  M  is  the  order  of  the  MA  part  and  A'  is  the 
order  of  the  AR  part.  In  contrast,  the  block  IV  and  ordinary7  least-squares  methods  re- 
quire computations  on  the  order  of  (M  +  A")2  . 

B.  PREVIOUS  MULTICHANNEL  METHODS 

Whittle  [Ref.  6:  pp.  129-130]  was  the  first  to  develop  a  multichannel  solution  for  the 
ARMA  modeling  problem.  He  sought  to  extend  the  recursive  Durbin  solution  for  esti- 
mating the  parameters  of  a  single  variable  autoregressive  process  to  a  multivariable 
autoregressive  process.  He  discovered  that  to  do  this  he  would  have  to  fit  the  data  to 
two  autoregressive  processes  simultaneously.  One  of  the  autoregressions  would  use 
present  data  samples  to  predict  the  value  of  the  data  one  time  step  in  the  future.  This 
is  called  forward  prediction.  The  second  autoregression  would  use  present  data  samples 
to  predict  the  value  of  the  data  at  the  previous  time  instant  and  is  referred  to  as  back- 
ward prediction.  Sometime  during  this  research,  Whittle  determined  that  if  the  input 
was  derived  from  a  MA  scheme,  making  the  process  ARMA,  then  the  solution  would 
remain  the  same  provided  the  correlations  of  the  input  used  in  the  parameter  estimation 
process  had  shifts  greater  than  the  MA  scheme.  Whittle's  use  of  the  two  separate  and 
simultaneous  autoregressions  to  model  an  ARMA  process  can  be  thought  of  as  a 
multichannel  modeling  approach. 

Further  work  in  the  area  of  multichannel  ARMA  modeling  was  conducted  by  Perry 
and  Parker  [Ref.  7:  pp.  509-510].  They  started  out  with  the  ARMA  problem  formulation 
discussed  in  Chapter  2.  Using  the  method  of  ordinary  least-squares  to  minimize  the 
mean  square  error,  they  found  the  solution  for  the  estimate  of  the  parameters  to  be  the 
Wiener  solution  given  by: 
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In  equation  (4.1)  Riy  is  a  matrix  of  autocorrelations  of  the  past  outputs.  RuV  is  a  matrix 
of  autocorrelations  of  the  inputs,  R>v  and  R^  are  crosscorrelations  of  the  input  and 
output  data,  b  is  a  vector  of  AR  parameters,  and  a'  is  a  vector  of  MA  parameters.  In 
addition,  r„  is  a  vector  of  autocorrelations  of  past  output  data  with  the  current  output 
and  r>v  is  a  vector  of  crosscorrelations  of  input  data  with  the  current  output.  By  as- 
suming the  first  MA  parameter,  a\  was  known,  they  were  able  to  treat  it  as  a  gain  and 
factor  it  out  of  all  the  other  MA  parameters.  This  allowed  them  [Ref.  7:  pp.  509-510] 
to  extract  the  (A'+  1)"  row  and  column  of  equation  (4.1)  and  rewrite  the  solution  in  the 
form: 
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In  equation  (4.2),  a'  has  been  rewritten  as  a  to  indicate  that  the  a'0  term  has  been  fac- 
tored out  of  all  of  the  MA  parameters.  Reasoning  that  equation  (4.2),  the  ARMA  sol- 
ution, was  a  generalization  of  the  AR  solution,  they  figured  that  it  must  have  a  recursive 
solution  consisting  of  some  combination  of  the  Levinson-Durbin  algorithm,  a  recursive 
solution  for  the  AR  problem.  They  then  determined  equation  (4.2)  was  in  a  form  similar 
to  Whittle's  formulation  of  the  problem.  So  they  reasoned  that  they  could  use  a  form 
of  Whittle's  solution  to  solve  the  ARMA  modeling  problem.  Like  Whittle,  their  solution 
consists  of  a  forward  and  a  backward  autoregression.  It  uses  two  coupled  lattice  filters 
to  process  the  input  and  output  data.  Off  diagonal  elements  of  the  lattice  coefficient 
matrices  specify  the  coupling  points  of  the  two  lattices. 

C.     ITERATIVE  APPROACH  TO  MULTICHANNEL  ARMA  MODELING 

This  thesis  proposes  another  solution  to  the  ARMA  modeling  problem  using  the 
multichannel  approach.  It  is  an  iterative  approach  with  no  direct  coupling  of  the  two 
channels.  However,  note  that  there  is  an  implicit  coupling  in  the  sense  that  the  ARMA 
system  output  samples  y{n)  are  a  function  of  the  present  and  past  input  samples  u(n), 
and  the  past  output  samples.  This  is  shown  in  equation  (2.1).  The  approach  proposed 
uses  two  separate  finite-impulse-response  (FIR)  filters  to  estimate  the  unknown  system. 
One  filter  estimates  the  AR  part  of  the  unknown  system.  The  second  estimates  the  MA 
part  of  the  unknown  system.   Figure  7  on  page  25  shows  the  structure  of  this  approach. 


23 


The  derivation  of  the  equations  for  this  approach  follows  the  method  of  ordinary  least- 
squares. 

A(z) 
From  Figure  7  on  pase  25,  the  value  of  the  sisnal  Y(z)  is  seen  to  be  ■        ■  Liz)  . 

B(z) 

When  this  signal  passes  through  C{z)  ,  if  C(z)  is  close  to  B(z),  the  resulting  signal  X(z)  is 
approximately  A(z)Lr(z).  Also  from  Figure  7  on  page  25,  the  value  of  Z(z)  is  seen  to  be 
A(z)U(z)  provided  D(z)  is  close  to  A(z).  The  difference  of  these  two  signals  forms  the 
error  which  we  seek  to  minimize  by  the  method  of  ordinary  least-squares.  In  minimizing 
the  error  we  seek  to  drive  D(z)  and  C(z)  as  close  as  possible  to  A(z)  and  B(z)  ,  respec- 
tively. 

Referring  to  Figure  7  on  page  25,  signals  z(n)  and  x(n)  are  defined  as  the  outputs 
from  two  FIR  filters  and  are  given  by  the  equations: 

z{n)  =  d0u{n)  +  d^u{n  -  1)  +  d2u{n  -  2)  H h  dMu{n  -  M)  =  uT{n)d  (4.3) 

x(n)  =  y(n)  +  cxy{n  -  1)  +  c^{n  -2)  +  -  +  c^y{n  -  N)  =  yr(«)c  (4.4) 

where  the  vectors  d  and  c  represent  the  systems  D(z)  and  C(z),  respectively.  The  vectors 
d,  u(/;).  c,  andy(//)  are  given  by  the  following  equations: 

d=ry0  dx  d2  ...  dMf  (4.5) 

u(«)  =  [u{n)    u{n-\)    u{n-2)     ...     u{n  -  M)f  (4.6) 

c  =  [l    q    c2     ...  ciV]r  (4.7) 

>'(")  =  bin)   y(n  -  1)   y{n  -  2)     ...y{n  -  N)f  (4.8) 

The  d  parameters  are  estimates  of  the  MA  portion  of  the  ARMA  process.  The  c  pa- 
rameters approximate  its  AR  portion.  The  vector  u(n)  is  the  input  data  vector  of  length 
M,  the  order  of  the  MA  part,  and  y(«)  is  the  output  data  vector  equal  of  length  Ar,  the 
order  of  the  AR  part. 

Forming  the  error  between  x  and  z  results  in: 

e(n)  =  z(n)  -  x(n)  =  uT(n)d  -  yT(n)c  (4.9) 

The  least-squares  performance  criterion  is: 
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Figure  7.      Multichannel  modeling  process 


L  L 


x(n)J 


(4.10) 


n=]  n=\ 


Substituting  for  e(n)  from  equation  (4.9),  we  can  write  equation  (4.10)  in  vector  form  as: 


J=  )  (urd-yrc)2 


n=\ 


(4.11) 


where  we  have  dropped  the  time  index,  n,  for  convenience.    Expanding  equation  (4.11) 
results  in: 


J  =   )  druurd  +  cryyrc  -  2druyrc 


(4.12) 


*=i 
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We  notice  that  the  performance  criterion  is  a  function  of  both  d  and  c.  Minimizing  the 
performance  criterion  by  differentiating  with  respect  to  the  vector  c  and  equating  the 
results  to  zero  yields: 


L  L 

M.  =  o  =  o  +  l)  (yyr)c  -  2)  (yur)d 


(4.13) 


n=\ 


Similarly,  differentiating  the  performance  criterion  with  respect  to  the  vector  d  and 
equating  the  result  to  zero  yields: 


cJ 
<5d 


0  =  0+2)  (uur)d  -  2  )  (uyr)c 


(4.14) 


n=\ 


n=\ 


Solving  equation  (4.13)  for  c  and  equation  (4.14)  for  d  results  in  two  equations  for 
estimating  the  AR  and  MA  parameters  of  the  unknown  system  given  by: 


C      =      Ryy     Ry^ 


(4.15) 


and 


d  =  RuX>C 


(4.16) 


where 


n=\ 


UO)    rJX) 


'uuiK) 
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(4.17) 


is  an  estimate  of  the  input  autocorrelation  matrix.    The  elements  of  Ruu  are  computed 
usine  the  unbiased  method  as  follows: 


L-l 


ruu{t)  =  -[zr[YJu{j)u{j-[)  for  /  =  0'  l>  2>  -  >M 


(4.18) 
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Matrices  Ryy,  R^,  and  R>a  appearing  in  equations  (4.15)  and  (4.16)  have  structures  iden- 
tical to  equation  (4.17),  where  Ryy  is  the  estimate  of  the  output  autocorrelation  matrix, 
and  RJ}  and  Ryu  are  estimates  of  the  crosscorrelation  matrices.  Note  that  Ry„  =  R^  .  The 
elements  of  these  matrices;   ryy,  ruy  and  ryu,  are  computed  as  follows: 

L-l 

^(O-T^ZvOW-O  for  /  =  0'  l'  2'  -  ,N  (419) 
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ruy 


(0  =y^  £"(/>(/ ~0  for  '=0.  1.  2,    ...  ,M  (4.20) 

7=0 


and 


L-l 

M/)  =  7TyZl'w^~/)  for  /=0'  l'  2'  -  'A'  (421) 

7=0 

Up  to  this  point,  following  the  standard  least-squares  procedure  has  resulted  in  two 
dependent  or  coupled  equations  to  solve  for  the  parameters  of  an  unknown  system 
modeled  as  an  ARMA  process.  How  best  to  solve  these  equations?  By  iteration.  The 
steps  of  the  iterative  process  are  to 

•  Calculate  the  correlation  matrices  and  vectors  from  the  available  data. 

•  Initialize  c  .    Exploit  the  fact  that  the  first  parameter  in  c  is,  c0  =  1. 

•  Solve  for  d  from  this  initial  c  . 

•  Solve  for  c  from  d  . 

•  Repeat  beginning  at  the  third  step. 

Here  is  a  summary  of  the  equations  in  proper  order  for  implementing  the  algorithm: 

•  Compute  Ruu  ,  R„  and  Ruy  from  equations  (4.17)  to  (4.21).   Note  that  R>u  =  R£. 

•  Initialization: 


For  k  =  0  toK 


where  k  is  the  iteration  number. 
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c(0)  =  [l    0     ...     0]7"  (4.22) 


d(*+1>  =  R;XvCW  (4-23) 

c(*+1)  =  R-1RVHd(*+1>  (4.24) 


This  is  a  very  simple  algorithm.  For  the  case  where  the  AR  and  MA  orders  are  equal, 
the  correlation  matrices  are  half  the  size  of  the  block  data  matrices  which  must  be  gen- 
erated and  inverted  in  the  IV  algorithm. 

1.     Testing  the  Multichannel  Iterative  Algorithm 

We  tested  the  algorithm  by  computer  simulation  using  second  and  third-order 
filters  as  unknown  systems.  Table  3  shows  pole  and  zero  locations  as  well  as  MA  and 
AR  parameters  for  the  test  filters.  Data  lengths  of  500  data  points  were  used  to  calculate 
the  autocorrelation  and  crosscorrelation  matrices.  Besides  the  reported  cases,  we  tested 
the  algorithm  on  several  other  second,  third  and  mixed-order  cases.  The  performance 
of  the  algorithm  in  all  cases  that  we  tested  was  fairlv  uniform. 


Table  3.     TEST  SYSTEMS  FOR  ITERATIVE  MULTICHANNEL  MODELING 
METHOD 


TEST 
FILTER 

location 
of  po:.es 

LOCATION 
OF  ZEROS 

AR 
PARAMETERS 

MA 
PARAMETERS 

T2 

0.445  +  J0.228 
0.445  -  J0.22S 

0.4  +  J1.273 

0.4  -jl.273 

1.0 

-0.S9 
0.25 

0.5 
-0.4 
0.89 

T2N 

0.445  +  J0.228 

0.445  -J0.228 

0.4  +  j0.8 
0.4  -J0.S 

1.0 
-0.89 
0.25 

1.0 
-0.80 

O.SO 

T3 

0.6605 
0.6647 +  J0.5020 

0.6647  -jo. 5020 

-1.0 
-1.0 
-1.0 

1.0 
-1.99 
1.572 

-0.45S3 

0.0154 
0.0462 
0.0462 
0.O154 

Results  of  the  tests  are  shown  in  graphical  form  in  Figure  8  on  page  29, 
Figure  9  on  page  30,  and  Figure  10  on  page  31.  Dashed  lines  indicate  the  true  values 
of  the  parameters.    Solid  lines  show  the  values  the  algorithm  calculated. 

Table  4  on  page  32  contains  the  algorithm's  best  estimates  of  the  parameters, 
along  with  the  number  of  iterations  required  to  converge  to  those  estimates.  It  also 
shows  the  absolute  and  percent  errors  from  the  true  parameters.  For  the  second-order 
test  cases,  convergence  to  the  true  parameter  values  occurred  within  20  iterations.  The 
third-order  test    case  took  14  iterations  to  converge  to  its  steady-state  values;  however, 

these  values  were  not  the  true  parameter  values. 
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Figure  8.      Second-order  test  case  T2.  (A)  MA  parameters.  (B)  AR  parameters. 
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Figure  9. 


Second-order  test  case  T2N.  (A)  MA  parameters.  (B)  AR  parameters. 
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Figure  10.      Third-order  test  case  T3.  (A)  MA  parameters.   (B)  AR  parameters. 
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Table  4.     PARAMETER  ESTIMATES  BY  THE  ITERATIVE  MULTICHANNEL 
ALGORITHM. 


TEST 

PARAMETER 

ABSOLUTE 

PERCENT 

ITERATIONS 

FILTER 

ESTIMATE 

ERROR 

ERROR 

T2 

0.500 

0.0 

0.0 

5 

-0.393 

-0.007 

1.75 

0.S90 

0.0 

0.0 

1.000 

0.0 

0.0 

-0.SS9 

+  0.001 

0.11 

0.247 

-0.003 

1.20 

T2X 

l.ooo 

0.0 

0.0 

20 

-0.794 

+  0.006 

0.75 

0.79S 

-0.002 

0.25 

1.000 

0.0 

0.0 

-0.SS6 

+  0.004 

0.45 

0.245 

-0.005 

2.00 

T3 

0.0153 

-0.0001 

0.65 

14 

0.0487 

+  0.0025 

5.41 

0.0590 

+  0.0128 

27.71 

0.0287 

+  0.0133 

86.36 

0.99 

-0.01 

1.0 

-1.79 

+  0.20 

10.05 

1.267 

-0.305 

19.40 

-0.3086 

+  0.1497 

32.66 

2.     Stopping  Parameter 

In  tests  of  third-order  systems,  the  parameter  estimates  swung  through  or  close 
to  the  true  coefficient  values  and  converged  to  values  somewhat  removed  from  the  true 
values.  We  developed  a  stopping  parameter  to  flag  the  iteration  when  the  estimates  were 
closest  to  the  true  values.  This  occurs  when  the  error  is  smallest.  Referring  to 
Figure  7  on  page  25.  if  D(r)  is  equal  to  A(z)  and  C(z)  is  equal  to  B(-),  x  and  z  will  both 
equal  A(r)U(z).  The  error  will  be  zero.  The  farther  removed  D(r)  and  C(z)  are  from 
A(r)  and  B(r),  respectively,  the  larger  the  error  becomes. 

The  stopping  parameter  is  calculated  from  the  difference  of  the  z  and  x  values 
at  every  iteration.  After  the  parameter  vectors  d  and  c  are  estimated  for  a  particular  it- 
eration, the  stopping  parameter  is  calculated  from: 


ek{n)  =  zk{n)  -  xk{n)  =  \Jxk  -  u[d/c 


(4.25) 
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where  k  is  the  iteration  number  and  d  .  u.  c  and  y  are  defined  in  equations  (4.5)  to  (4. Si. 
The  parameter  vectors  d  and  c  represent  the  systems  D(r)  and  C{z),  respectively. 

Figure  11  on  page  34  and  Figure  12  on  page  35  graph  the  stopping  parameter 
(dotted  line)  along  with  the  estimated  and  true  values  of  the  parameters.  They  show  that 
when  the  stopping  parameter  is  smallest,  the  parameters  are  closest  to  their  true  values. 
Table  5  shows  the  improvement  in  the  estimates  of  the  parameters  resulting  from 
choosing  those  values  when  the  stopping  parameter  is  smallest. 


Table  5.     PARAMETER    ESTIMATES   WHEN   STOPPING   PARAMETER    IS 
SMALLEST. 


TEST 

PARAMETER 

ABSOLUTE 

PERCENT 

ITERATIONS 

FILTER 

ESTIMATE 

ERROR 

ERROR 

T3 

0.0156 

+  0.0002 

1.30 

10 

0.0478 

+  0.0016 

3.46 

0.0536 

+  0.0074 

16.02 

0.0196 

+  0.0042 

27.27 

0.97 

-0.03 

3.0 

-1.84 

+  0.15 

7.54 

1.375 

-0.197 

12.53 

-0.3672 

+  0.0911 

19.88 

The  stopping  parameter  can  be  used  in  a  real  modeling  situation  because  it 
comes  from  the  data  and  the  estimates  of  the  parameters.  Another  measure  of  how  well 
the  estimates  of  the  parameters  fit  the  actual  system  is  the  norm  of  the  coefficient  error. 
This  cannot  be  used  in  a  real  modeling  situation  however,  because  the  values  of  the  true 
parameters  are  not  known.  We  calculated  it  for  the  test  cases  as  a  check  on  the  appro- 
priateness of  using  the  stopping  parameter.  Figure  13  on  page  36  and  Figure  14  on 
page  37  graph  the  stopping  parameter  (dotted  line)  and  the  norm  of  the  coefficient  error 
for  test  cases  T2  and  T3.  On  both  graphs  the  two  curves  correspond  well.  Both  reach 
their  minimum  value  at  the  same  point,  the  point  where  the  estimates  of  the  parameters 
are  closest  to  their  true  values. 

3.     Linear-prediction  of  the  Denominator  Coefficients 

The  iterative  algorithm  detailed  in  equations  (4.22)  to  (4.24)  starts  by  initializing 
the  AR  parameter  estimates  to: 


c(0)  =  [l    0     ...     0]r 


(4.26) 
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Figure  11.      Stopping  parameter  example  for  test  case  T2. 
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Figure  12.      Stopping  parameter  example  for  test  case  T3. 
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Figure  13.      Stopping  parameter  and  norm  of  the  coefficient  error  for  test  case  T2. 

where  c0  is  known  to  be  1  from  the  Z  transform  of  the  ARMA  difference  equation, 
equation  (2.1).   The  Z  transform  is  given  by: 

¥(-)[{  +  cjz"1  +  cV"2  +  -]=  l'(--)M)  +  4-"1  +  ^'"2  +  -  1  (4-7) 


-i 


Fir)        ^j  +  ^/,2      +d2:   *  + 


t/(z] 


1  +  C,Z        -I-  C-,1        + 


(4.28) 


The  initial  estimates  for  the  other  AR  parameters  are  zero.  This  can  be  far  from 
their  actual  values.  A  closer  estimate  of  the  other  AR  parameters  should  result  in 
quicker  convergence  for  all  parameters.  A  closer  estimate  of  the  AR  parameters  can  be 
obtained  by  using  linear-prediction  techniques.  Figure  15  on  page  39  shows  the  ap- 
proach used.  In  Figure  15  on  page  39,  y(n)  is  the  output  from  the  unknown  system. 
The  system  C'(z),  which  is  represented  by  the  vector  c'  ,  is  the  linear-prediction  filter  used 
to  estimate  >•(«).  It  uses  the  previous  n  —  S  samples  of  the  output  to  generate  a  current 
estimate.   This  is  given  by  the  equation: 
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Figure   14.      Stopping  parameter  and  norm  of  the  coefficient  error  for  test  case  T3. 


y{n)  =  yr(«-l)c'  (4.29) 

where  c'  is  a  vector  of  the  tap  weights  of  the  autoregressive  process  given  by: 

c'  =  [c'0    c\     ...c'„_v]r  (4.30) 

and  y(/7  —  1)  is  a  vector  of  the  past  output  values  given  by: 

y(«  -  1)  =  \y{n  -  1)   y{n  -  2)     ...     y{n  -  N)f  (4.31 ) 

Following  least-squares  techniques,  we  form  the  error  between  the  estimate  and  the  ac- 
tual value  of  the  output.  The  sum  of  the  squares  of  the  errors  becomes  the  performance 
criterion.  This  is  differentiated  with  respect  to  the  tap  weights  and  set  equal  to  zero. 
Solving  this  for  the  tap  weights  results  in  the  equation: 


C'     -     Ryy    Tyy 


(4.32) 
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This  is  the  standard  Weiner  filter  solution  [Ref.  8:  p.  32].  It  tells  us  that  the  best  estimate 
of  the  AR  parameters  can  be  found  from  the  correlations  of  the  output  data.  The  matrix 
R(>  is  the  autocorrelation  matrix  of  the  past  outputs  and  ryy  is  autocorrelation  vector  of 
the  past  outputs  with  the  current  output.  In  all  cases  tested,  we  did  not  achieve  any 
significant  improvement  in  the  accuracy  of  the  estimates  of  the  parameters,  or  in  the 
speed  of  convergence,  using  the  straight  linear  prediction  of  equation  (4.32). 

A  modification  to  this  approach,  which  we  refer  to  as  modified  linear-prediction, 
uses  correlation  lags  beginning  on  the  order  of  the  MA  portion  of  the  ARMA  process. 
For  example,  correlations  for  calculating  Ryy  for  a  third-order  system  would  start  at  a  lag 
of  three  and  increase  to  a  lag  of  five.  Correlations  for  calculating  ryy  would  start  at  a  lag 
of  four  and  increase  to  a  lag  of  six.  This  ensures  that  the  effect  of  the  MA  part  of  the 
unknown  system  is  removed  from  the  linear-prediction  of  the  AR  part.  This  modified 
method  of  linear-prediction  is  given  by  the  equation: 


n-X 


ryy{q)  ryy{q-\) 

ryy{q+\)  ryy(q) 


rvy{q  +  p-  1)    ryy{q  +  p-2) 
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ryy{q  - 

-p  + 

1) 

ryy(q  - 

-P  + 

2) 

ryy(q) 

ryy(q+  1) 
ryy{q  +  2) 


'Vy(<7  +  P) 


(4.33) 


where  q  is  the  order  of  the  MA  portion  and  p  is  the  order  of  the  AR  portion.    [Ref.  2: 
p.  182] 

Figure  16  on  page  40  shows  the  results  of  using  modified  linear-prediction  with 
third-order  test  case  T3.  When  comparing  this  graph  to  the  estimates  obtained  without 
linear-prediction,  shown  in  Figure  12  on  page  35,  note  that  the  vertical  axes  have  dif- 
ferent scales.  Table  6  on  page  39  lists  the  values  of  the  estimates  at  iteration  10  and 
compares  them  with  the  true  values  via  the  absolute  and  percent  errors.  A  comparison 
of  Table  6  on  page  39  with  Table  5  on  page  33,  the  best  estimates  without  the  use  of 
modified  linear-prediction,  shows  that  modified  linear  prediction  has  significantly  in- 
creased the  accuracy  of  the  AR  estimates  at  the  tenth  iteration.  The  accuracy  of  the 
MA  estimates  remains  approximately  the  same.  The  tenth  iteration  was  chosen  as  the 
point  to  select  the  parameter  values  because  in  both  cases  this  was  the  iteration  where 
the  stopping  parameter  had  the  smallest  value. 
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Figure   15.      Linear  prediction  block  diagram 


Table  6.     PARAMETER  ESTIMATES  USING  MODIFIED 

LINEAR-PREDICTION  FOR  INITIAL  ESTIMATE  OF  AR  PARAME- 
TERS. 


TEST 

PARAMETER 

ABSOLUTE 

PERCENT 

ITERATIONS 

FILTER 

ESTIMATE 

ERROR 

ERROR 

T3 

0.015(3 

+  0.0002 

1.30 

10 

0.0485 

+  0.0023 

4.98 

O.D512 

+  0.0050 

10.80 

0.0101 

+  0.0053 

34.42 

LOO 

0.0 

0.0 

-1.98 

+  0.01 

0.50 

1.553 

-0.019 

1.21 

-0.4458 

+  0.0125 

2.73 

D.     FORMULATION  OF  THE  SEQUENTIAL  MULTICHANNEL  APPROACH 

To  decrease  the  computational  intensity  of  updating  the  estimates  of  the  AR  and 
MA  parameters  due  to  new  data,  an  algorithm  to  sequentially  process  the  data  has  been 
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Figure  16.      Parameter  estimates  for  test  case  T3  using  modified  linear-prediction  for 
initial  estimate  of  AR  parameters. 
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developed.  Development  begins  with  the  performance  criterion  seen  previously  for  the 
block  data  case  in  equation  (4.10).  From  that  starting  point,  equations  are  developed 
which  relate  new  estimates  of  the  parameters  to  the  previous  estimates  and  the  new  data. 
Separate  but  similar  equations  are  developed  for  the  MA  and  the  AR  coefficients.  Due 
to  the  similar  nature  of  the  development  of  these  equations,  only  the  development  of  the 
equations  for  the  AR  coefficients  is  presented  here.  However,  the  final  results  for  both 
the  AR  and  MA  coefficients  are  given. 

The  performance  criterion  can  be  written  as: 


n  n 


(4.34) 


(=0  (=0 


Expanding  this  results  in: 


J  =  V z[z{  -  2z[y[c  +  Jyxfc  (4.35) 

where  c  and  y  are  defined  in  equations  (4.7)  and  (4.8)  and  z  is  the  scalar  signal  at  the 
output  of  D.  Differentiating  the  performance  criterion  with  respect  to  c  and  setting  the 
result  equal  to  zero  yields: 

n  n 

4^  =  o  =  £>f)c-;>><  (4.36) 

j=0  i"=0 


Equation  (4.36)  can  also  be  written  as 

n  n 

/=0  1=0 

Solving  for  the  AR  parameter  vector  results  in: 

n  n 


i=0  (=0 
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Because  c  is  an  estimate  based  on  data  available  through  time  n  we  signify  this  by  in- 
troducing the  index  n  on  c  to  yield: 


/=0  /=0 


The  first  step  in  formulating  the  sequential  algorithm  is  to  develop  an  update 
equation  for  the  estimate  of  the  AR  parameters.  Applying  the  method  presented  in 
Graupe  [Ref.  9:  p.  124],  we  first  define  a  new  matrix  P~l  as 


ri 

^-£yjr  (4.40) 


1=0 

This  is  a  matrix  of  the  output  data  of  the  unknown  system.   At  the  previous  time  n  —  1 
this  matrix  is  written  as: 

M-l 

P*li  =  £y,yr  (4.4i) 

i=0 

By  substituting  equation  (4.40)  into  equation  (4.39)  the  estimate  of  the  AR  parameters 
can  be  rewritten  as: 

n 

cn  =  P^W  (4-42) 

/=o 

The  right  side  of  equation  (4.42)  needs  to  be  converted  into  an  expression  containing  the 
previous  estimate  of  the  parameters  plus  a  correction  term.  It  needs  the  past  value  of 
c„  which  is  c„_,.    From  equation  (4.39)  we  can  write: 

n-]  H-l 

c„_,  =  (Xy/y/rr!Sz^  (443) 

i=C  /=0 


This  can  be  rewritten  as: 


n-\  n-\ 


X^- = (Z>'^c-i  (444^ 


j'=0  /=0 
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Premultiplying  equation  (4.42)  by  P;1  and  then  separating  the  last  term  from  the  sum- 
mation results  in: 


n-\ 


pA„= 2/^+2^ '  (4-45) 


(=0 


.=0 


Substituting  for  £z,y,  in  equation  (4.45)  from  its  equivalent  expression  in  equation  (4.44) 
vields: 


n-\ 


p«1c=(2Jy/y/r)cn-i+^y«  (4-46) 

By  adding  y„y^cn_,  to  the  right  side  of  equation  (4.46)  and  grouping  it  with  the  summa- 
tion, the  summation  will  range  from  /  =  0  to  n  .  In  order  not  to  affect  the  value  on  the 
right  side  of  equation  (4.46),  y„ynrc„-i  must  also  be  subtracted  from  the  right-hand  side, 
which  yields: 

p^c  =  (2jfyf)c„-i  +  z„y„  +  y«y^M-i  -  y^^-i  (4.47) 

Combining  ynyjcn_,  with  the  summation  as  describe  above  results  in: 

n 

P;l,c  =  (2]y/yr)cn_,  +  inyn  -  ynyjew_,  (4.48) 


(=0 


Replacing  £yfyf  with  its  equivalent  expression  from  equation  (4.40)  yields: 


1=1 


P«  lc«  =  P*  lc„_,  +  y„(z„  -  y„Vi)  (4-49) 

Premultiplying  by  P„  results  in  the  following  equation  for  the  update  of  the  estimate  of 
the  AR  parameters: 

c„  =  c„_,  +  ?nyn(zn  -  yjcn_,)  (4.50) 

This  is  the  desired  result.    It  relates  the  estimate  of  the  parameters  at  time  A'  to  the  es- 
timate at  the  previous  time,  A'—  1,  plus  the  new  output  data  vector,  y„,  and  the  error  at 
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time  X.  Error  is  represented  by  the  z„  —  y;[cn_,  term.   The  corresponding  equation  for  the 
MA  parameters  is: 

d„  =  d„_!  +  Q/2u„(x„  -  u„V,)  (4.51) 

In  equation  (4.51),  Q  is  a  matrix  of  the  input  data  of  the  unknown  system  given  by: 

n 

q~'  =  £u.uf  (452) 

;=o 

Finally,  we  need  a  sequential  update  for  P„  and  Q„  to  complete  the  sequential  algo- 
rithm.  This  is  accomplished  by  using  a  form  of  the  matrix  inversion  lemma. 

By  pulling  the  last  term  out  of  the  summation,  the  definition  of  P~l  given  in  equation 
(4.40)  can  be  rewritten  as: 

n-\ 

Pn]=Y,y^+ynJn  (4-53) 

(=0 

n-l 

Substituting  for  Zy,y,r  its  equivalent  expression  from  equation  (4.41)  results  in 

1=0 

p;1  =  p£i  +  ynyn  (4-54) 

Inverting  both  sides  of  the  equation  results  in: 

Pn  =  (P-1,  +  y^y7)"'  (4.55) 

Let  A  =  P-J,,  B  =  y„,  C  =  1  ,  and  D  =  yj.   Then,  by  the  matrix  inversion  lemma: 

P„  =  A-1  -  A-1B(C_1  +  DA~1B)"1DA~1  (4.56) 

Substituting  the  appropriate  expressions  into  equation  (4.56)  results  in: 

p„  =  (P^-i)"1  -  (P^i)_1y,[i  +  yr(P;l1r'yJ"1yr(P;l,r1  (4.57) 


This  reduces  to: 


p„  =  p„_t  -  Pn_iyn[l  +  y„rPfl.1yn]-1y„rP„-,  (4-58) 


Using  this  same  procedure,  the  update  for  Q„  is: 
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Q«  -  Q„-i  -  Qn-iUnU  +  u„rQ„-iU„]   ]^Qn_,  (4.59) 

A  reduction  in  the  computational  intensity  has  been  achieved  by  reducing  the  matrix 
inversion  of  equation  (4.55)  to  inversion  of  a  scalar  in  equation  (4.57).  Inversion  of 
these  scalars  is  much  simpler  than  inversion  of  the  Ryy  and  Ruu  matrices  of  the  block 
processing  case. 

The  sequential  multichannel  algorithm  is  summarized  below: 

•  The  parameter  update  equations: 

c„  =  c„_j  +  Pnyn{zn  -  yjcn-i)  (4.60) 

d„  =  d„_,  +  Q„u„(j:n  -  uX-i)  (4.61) 

•  The  update  equations  for  the  P  and  Q  matrices: 

p«  =  Pn-i  -  P,-iy.[i  +  yjp^iyj"!yfo-i  (4-62) 

Q„  =  Q,?_,  -  Q„_!Un[l  +  ujQ^iUj^nJOi-i  (4.63) 

The  reduction  in  computational  intensity  comes  with  a  trade-off.  Now  the  algo- 
rithm is  more  complex  to  use.  Updates  are  required  for  P  and  Q  as  well  as  c  and  d  where 
before,  in  the  block  multichannel  algorithm,  updates  were  only  required  for  c  and  d.  But. 
as  in  the  sequential  IV  case,  an  added  advantage  of  the  sequential  multichannel  algo- 
rithm is  it  allows  updates  of  the  estimates  of  the  parameters  based  upon  new  data. 
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V.     SUMMARY 

In  this  thesis  we  set  out  to  develop  two  algorithms  for  modeling  unknown  systems 
as  ARMA  processes.  These  are  the  IV  method  of  system  identification  presented  in 
Chapter  3,  which  is  a  modification  of  the  method  of  ordinary  least-squares,  and  the  it- 
erative multichannel  method  presented  in  Chapter  4. 

The  IV  method  is  not  a  new  concept  in  either  its  block  or  sequential  processing 
forms.  However,  our  derivation  of  the  sequential  algorithm  was  done  independently  of 
other  IV  sequential  algorithms.  We  chose  the  IV  method  because  it  reportedly  has  good 
noise  handling  capabilities  and  yields  consistent  and  unbiased  estimates  of  the  unknown 
system's  parameters.  It  also  remains  as  easy  to  use  as  the  method  of  ordinary  least- 
squares.  We  also  wanted  to  gain  familiarity  with  it  because  it  was  a  possible  candidate 
for  incorporation  into  the  multichannel  method. 

Operating  alone,  the  IV  method  produces  accurate  estimates  of  the  unknown  sys- 
tem's parameters.  Convergence  was  within  20  iterations  for  several  second-order  sys- 
tems that  we  tested.  Convergence  slows  down  as  the  system  order  increases.  However, 
the  results  do  converge  to  the  actual  system  parameters  given  sufficient  processing  time. 
The  performance  of  the  IV  method  is  similar  to  the  performance  of  the  method  of  ordi- 
nary least-squares. 

The  proposed  iterative  multichannel  algorithm  is  new  in  both  its  block  and  sequen- 
tial processing  forms  to  the  best  of  our  knowledge.  It  is  very  simple  to  use  in  the  block 
form.  It  achieves  accurate  results  for  second-order  systems  but  worse  results  for  third- 
order  systems  with  block  correlation  elements  calculated  based  on  only  500  data  points. 
Implementing  the  stopping  parameter  increases  the  accuracy  when  the  parameter  esti- 
mates converge  but  not  to  the  true  parameters.  Due  to  its  ability  to  separately  process 
the  input  and  output  data  from  the  unknown  system,  correlation  matrices  in  the  multi- 
channel block  processing  case  are  half  the  size  of  correlation  matrices  required  for  the 
single  channel  block  processing  case.  This  feature  reduces  the  computational  intensity 
over  what  is  required  for  the  conventional  least-squares  processing  case.  The  number 
of  iterations  required  for  convergence  seems  to  be  independent  of  the  order  of  the  sys- 
tem. However,  the  accuracy  of  the  estimates  suffer  as  the  order  of  the  system  increases. 
Using  linear  prediction  to  estimate  the  initial  values  of  the  AR  parameters  did  not  speed 
up  convergence  or  increase  the  accuracy  of  the  parameter  estimates.    However,  using 
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modified  linear  prediction  significantly  increased  the  accuracy  of  the  AR  parameter  es- 
timates, although  it  had  no  effect  on  the  MA  parameter  estimates. 

We  formulated  the  multichannel  sequential  algorithm.  This  allows  the  estimates  of 
the  parameters  to  be  updated  as  new  data  becomes  available.  But  we  have  not  tested 
this  algorithm.  It  needs  checking  using  a  variety  of  second  and  third-order  test  cases  to 
verify  that  it  works.  During  testing,  guidelines  need  to  be  developed  for  the  best  meth- 
ods to  initialize  the  P  and  Q  matrices  to  achieve  the  quickest  convergence  and  the  most 
accurate  parameter  estimates. 

As  mentioned  above,  one  of  the  reasons  for  investigating  the  IV  method  of  system 
identification  was  for  possible  inclusion  into  the  multichannel  algorithm,  the  hope  being 
that  the  favorable  noise  performance  of  the  IV  method  would  improve  the  performance 
of  the  multichannel  method.   This  is  another  area  that  remains  unexplored. 

The  block  multichannel  and  IV  methods  achieved  similar  results  for  second-order 
test  cases.  Convergence  to  the  actual  system  parameters  came  within  20  iterations  for 
both  algorithms.  However,  for  third-order  systems,  convergence  was  much  quicker  with 
the  multichannel  block  algorithm  than  with  the  sequential  IV  algorithm.  But  the  pa- 
rameter estimates  by  the  IV  method  were  more  accurate  than  by  the  multichannel  block 
method.  A  combination  of  the  two  algorithms  has  the  potential  for  incorporating  their 
unique  advantages  into  a  better  overall  parameter  estimation  method. 

Areas  for  further  research  are  listed  below: 

•  Verify  that  the  multichannel  sequential  algorithm  developed  in  Chapter  4  works  as 
a  means  of  modeling  an  unknown  system  as  an  ARMA  process. 

•  Investigate  the  possibility  of  incorporating  the  IV  method  into  the  multichannel 
sequential  algorithm. 

•  Analyze  why  initializing  the  AR  parameters  to  the  values  calculated  by  linear  pre- 
diction improves  the  speed  of  convergence  of  the  AR  parameters  in  the  multi- 
channel block  algorithm  but  does  not  improve  the  convergence  of  the  MA 
parameters.  Identify  a  method  for  obtaining  an  initial  estimate  of  the  MA  param- 
eters to  improve  their  speed  of  convergence  and  accuracy. 

•  Investigate  the  effects  of  increasing  the  number  of  data  points  used  to  calculate  the 
correlation  matrices  for  the  multichannel  block  algorithm  on  the  accuracy  of  the 
parameter  estimates  and  their  speed  of  convergence. 

•  Investigate  the  performance  of  the  IV  method  with  noise  present.  Compare  this 
performance  with  the  performance  of  the  method  of  ordinary  least-squares  with 
noise  present. 
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APPENDIX 

A.     INSTRUMENTAL  VARIABLE  ALGORITHM 
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*  WSIZE  ORDER  OF  AR  PART  OF  THE  AUXILIARY  MODEL 

*  YSIZE  ORDER  OF  THE  AR  PART  OF  THE  TEST  SYSTEM 

*  USIZE  ORDER  OF  THE  MA  PART  OF  THE  TEST  SYSTEM  AND  AUXILIARY 

*  MODEL 

*  VARIABLES  THAT  END  IN  R  CONTAIN  THE  ROW  SIZE  OF  THEIR  RESPECTIVE 

*  MATRICES.   VARIABLES  THAT  END  IN  C  CONTAIN  THE  COLUMN  SIZE  OF 

*  THEIR  RESPECTIVE  MATRICES. 

****    VARIABLE  DECLARATIONS    **** 

COMMON  /D/  RIVCOF(0:  1000,10) 

REAL  Z(10,10),Z  TPO(10,10),X(10,10),X  TPO(10,10) 
REAL  U(1),Y(10,10),W(10,10) 
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REAL  QMAT(10,10),Q  TEMP( 10, 10) ,QTEMP2( 10 , 10) ,QTEMP3( 10, 10) 
REAL  IV( 10, 10) , IVTEMPC 10, 10) ,COR( 10, 10) ,COEF( 10, 10) 
REAL  SCALAR, SCALR2 

DOUBLE  PRECISION  SEED 

INTEGER  I,J,K,VSIZE,YSIZE,USIZE 

INTEGER  ZMR , ZMC , ZTR , ZTC , XMR , XMC , XTR , XTC 

INTEGER  UMR,UMC,YMR,YMC,WMR,WMC 

INTEGER  QMR , QMC , QTR , QTC , QT2R , QT2C , QT3R , QT3C 

INTEGER  I VR , I VC , I VTR , I VTC , CMR , CMC , COEFMR , COEFMC 

INTEGER  IVCR.IVCC 

LOGICAL  EOF 

READ(4,*,END=22)  YS I ZE, US I ZE , COEFMR, COEFMC, ITERA 
CALL  RDMAT(COEF, COEFMR, COEFMC) 

*  INITIALIZE  VARIABLES 

EOF  =  .FALSE. 
XTR  =  COEFMC 
XTC  =  COEFMR 
ZTR  =  COEFMC 
ZTC  =  COEFMR 
IVR  =  COEFMR 
IVC  =  COEFMC 
QMR  =  COEFMR 
QMC  =  COEFMR 
SEED  =  888. 0D1 
IVCR  =  0 
IVCC  =  COEFMR 
WSIZE  =  YSIZE 

*  ZERO  OUT  THE  IV  PARAMETER  VECTOR,  THE  AUXILIARY  MODEL  DATA  VECTOR 

*  AND  THE  TEST  SYSTEM  DATA  VECTOR. 

CALL  INIT(IV,IVR,IVC,0.0) 
CALL  INIT(Z  TPO, ZTR, ZTC, 0.0) 
CALL  INIT(X  TPO, XTR, XTC, 0.0) 

*  INITIALIZE  THE  QMAT  AS  A  DIAGONAL  MATRIX  WHOSE  DIAGONAL  ELEMENTS 

*  EQUAL  100 

CALL  IN1TD(QMAT, QMR, QMC, 100.0) 

*  GET  VALUE  FOR  U(0).   U  IS  A  GAUSSIAN  RANDOM  VARIABLE. 

CALL  GGNML  (SEED,1,U) 

*  SHIFT  U(0)  INTO  X  &  Z  VECTORS  TO  CREATE  X(0)  S.   Z(0) 

Y(l,l)  =  0.  0 

CALL  SHIFT(X  TPO, XTC, YSIZE, USIZE, Y( 1 , 1) ,U( 1)) 

W(l,l)  =0.0 

CALL  SHIFT(Z  TPO, ZTC, WSIZE , USIZE, W( 1 , 1) ,U( 1) ) 
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*  CALCULATE  Y(0)   =  X  TPO(O)   *  COEFFICIENT  VECTOR 

CALL  MULTI(X  TPO,XTR,XTC,COEF,COEFMR,COEFMC,Y, 1 , 1) 

*  CALCULATE  W(0)   =  Z  TPO(O)   *   IV  VECTOR 

CALL  MULTI(Z  TPO,ZTR,ZTC, IV, IVR, IVC ,W, 1, 1) 
DO   90   I   =  0,ITERA 

CALL  GGNML(SEED,1,U) 

CALL  TPOSE( IV , IVR , IVC , IVTEMP , IVTR , IVTC ) 

*  SAVE  THE   IV  PARAMETERS 


91 


DO   91   J  =   1,IVTC 

RIVCOF(I,J)   =   IVTEMP(1,J) 
CONTINUE 
CALL  PRMAT( IVTEMP, IVTR, IVTC) 


*  SHIFT  Y(M)   AND  U(M+1)    INTO  X  TPO(M)   TO  GET  X  TP0(M+1) 

CALL  SHIFT(X  TPO,XTC,YSIZE,USIZE, Y( 1 , 1) ,U( 1)) 

*  SHIFT  W(M)   AND     U(M+1)    INTO   Z  TPO(M)   TO  GET  Z  TP0(M+1) 

CALL   SHIFT(Z   TPO,ZTC ,WSIZE,USIZE ,W( 1 , 1) ,U( 1) ) 

*  CALCULATE  Y(M+1)   AND  Z(M+1) 

CALL  MULTI(X  TPO,XTR,XTC,COEF,COEFMR,COEFMC,Y, 1 , 1) 
CALL  MULTI(Z  TPO,ZTR,ZTC,IV, IVR, IVC,W, 1 , 1) 

*  CALCULATE  THE  NEW  Q  MATRIX 

CALL  MULTI(X  TPO , XTR , XTC , QMAT , QMR , QMC , Q  TEMP,QTR,QTC) 
CALL  TPOSECZ  TPO,ZTR, ZTC,Z ,ZMR, ZMC) 

CALL  CORE (TMAT, QMR, QMC, Z,ZMR, ZMC, X  TPO, XTR, XTC, COR, CMR, CMC) 
CALL  MULTI(COR,CMR,CMC,Q  TEMP,QTR,QTC ,QTEMP2,QT2R,QT2C) 
CALL  SUBTRC ( QMAT , QMR , QMC , QTEMP2 , QT2R , QT2C , QMAT , QMR , QMC ) 

*  CALCULATE  THE  NEW   IV  VECTOR 

CALL  MULTI(X  TPO, XTR, XTC , IV, IVR, IVC , IVTEMP, IVTR, IVTC) 

SCALR2  =  Y(l,l)    -    IVTEMP(1,1) 

CALL   SMULTI ( SCALR2 , COR , CMR , CMC ) 

CALL  ADD(IV,IVR,IVC,COR,CMR,CMC,IV,IVR,IVC) 


90 


CONTINUE 


*  PLOT  THE   IV  PARAMETERS  VS  THE   ITERATION  NUMBER 
CALL  PLOT2(ITERA,USIZE,YSIZE,COEF) 

22  STOP 

END 

ft  ftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftftft 

SUBROUTINE   C0RE(MAT1 , I1R, I1C,MAT2, I2R, I2C,MAT3, I3R,I3C,RMAT,IRR, 
+IRC) 
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REAL  MAT1(10,10),MAT2(10,10),MAT3(10,10),RMAT(10,10) 
REAL  Q  TEMP(10,10),QTEMP2(10,10) 
INTEGER  IRR , IRC , QTR , QTC , QT2R , QT2C 

*  CALCULATE  THE  CORE:  Qf M)Z(M+1) ° 1+X' (M+1)Q(M)Z(M+1)  **-l 

*  MAT1  IS  THE  Q  MATRIX,  MAT2  IS  THE  Z  VECTOR,  AND  MAT3  IS  THE 

*  X  VECTOR. 

CALL  MULTI(MAT1,I1R,I1C,MAT2,I2R,I2C,Q  TEMP , QTR , QTC ) 

CALL  MULTI(MAT3,I3R,I3C,Q  TEMP, QTR, QTC, QTEMP2,QT2R,QT2C) 

SCALAR  =  1/(1  +  QTEMP2(1,1)) 

CALL  SMULTK  SCALAR, Q  TEMP , QTR , QTC ) 

CALL  ADD(Q  TEMP , QTR , QTC , Q  TEMP , QTR, QTC ,RMAT, IRR, IRC) 

CALL  SMULTI(0.5,RMAT,IRR,IRC) 

RETURN 

END 

SUBROUTINE  PLOT2( ITER A , ICURVN , ICURVD , MAT1 ) 

COMMON  /D/  RIVCOF(0: 1000,10) 

REAL  X(0:  1000), Y(0:  1000)  , MATH  10, 10)  , MAX, MIN 

INTEGER  I,J,ITERA,ICURVN,ICURVD,ISTP 

CALL  LIMITS( ICURVN, ICURVD, NMAX,NMIN,NSTP, 
+DMAX,DMIN,DSTP,ITERA) 

*  GENERATE  THE  ITERATION  NUMBER 

DO  90  I  =  0,ITERA 

X(I)  =  I 

Y(I)  =  0.  0 
90    CONTINUE 

*  CALCULATE  X  AXIS  LABELING  INTERVAL 

ISTP  =  ITERA/10 

*  SET  UP  THE  PLOT 

CALL  SHERPA( ' IVGRAF   ',*A',3) 

CALL  RESET('ALL') 

CALL  PAGE(8.5  11.0) 

CALL  HEIGHT(0. 14) 

CALL  HWROT('AUTO') 

CALL  XINTAX 

CALL  YAXANG(O) 

CALL  PHYS0R(  1.5,6.  0) 

CALL  AREA2D(5.0,3.5) 

CALL  COMPLX 

CALL  XNAME( ' ITERATIONS$ ' , 100) 

CALL  YNAME(' COEFFICIENT  VALUES$ ' , 100) 

CALL  MESSAG('(A)$' ,100,2.4,-0.8) 

CALL  THKFRM(0.03) 

CALL  FRAME 


IVA01970 

IVA02000 

IVA02010 

IVA02020 

IVA02030 

IVA02040 

IVA02060 

IVA02070 

IVA02090 

IVA02100 

IVA02110 

IVA02120 

IVA02130 

IVA02140 

IVA02150 

IVA02170 

IVA02180 

IVA02190 

IVA03720 

IVA03730 

IVA03740 

IVA03750 

IVA03760 

IVA03780 

IVA03800 

IVA03810 

IVA03820 

IVA03830 

IVA03850 

IVA03860 

IVA03870 

IVA03880 

IVA03890 

IVA03900 

IVA03910 

IVA03920 

IVA03930 

IVA03950 

IVA03960 

IVA03970 

IVA03980 

IVA04010 

IVA04030 

IVA04040 

IVA04050 

IVA04060 

IVA04070 

IVA04080 

IVA04090 

IVA04100 

IVA04110 

IVA04140 

IVA04150 

IVA04160 

IVA04170 

IVA04180 
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CALL  GRAF(0,ISTP,ITERA,NMIN,NSTP,NMAX) 

*  PLOT  THE  NUMERATOR  VALUES 

DO  93  J  =  ICURVD  +  1,ICURVN  +  ICURVD 
DO  94  I  =  0,ITERA 

Y(I)  =  RIVCOF(I,J) 

94  CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
93    CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  COEFS'  TRUE  VALUES 

CALL  DASH 

*  PLOT  NUMERATOR  COEFS '  TRUE  VALUES 

DO  95  K  =  ICURVD  +  1, ICURVD  +  ICURVN 
DO  96  J  =  0,ITERA 
Y(J)  =  MAT1(K,1) 
96       CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 

95  CONTINUE 
CALL  ENDGR(O) 

*  SET  UP  SECOND  PLOT  FOR  DENOMINATOR  PARAMETERS 

CALL  RESET('DASH') 

CALL  PHYSOR(1.5,1.5) 

CALL  AREA2D(5.0,3.5) 

CALL  COMPLX 

CALL  XNAME( ' ITERATIONS$ ' , 100) 

CALL  YNAMEC PARAMETER  VALUES$ ' , 100) 

CALL  MESSAG('(B)$' ,100,2.4,-0.8) 

CALL  THKFRM(0.03) 

CALL  FRAME 

CALL  GRAF(0,ISTP,ITERA,DMIN,DSTP,DMAX) 

*  PLOT  THE  DENOMINATOR  VALUES 

DO  91  J  =  1, ICURVD 
DO  92  I  =  0,ITERA 

Y(I)  =  -RIVCOF(I,J) 
92       CONTINUE 

CALL  CURVE (X,Y,ITERA,0) 
91    CONTINUE 

*  PLOT  DENOMINATOR  COEFS '  TRUE  VALUES 

CALL  DASH 

DO  99  K  =  1, ICURVD 
DO  100  J  =  0,ITERA 
Y(J)  =  -MAT1(K,1) 
100      CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
99    CONTINUE 
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CALL  ENDPLCO) 
CALL  DONEPL 
RETURN 
END 

SUBROUTINE  SUBTRC(MAT1 , I 1R, I1C,MAT2 , I2R, I2C ,RMAT, IRR, IRC) 

*  **5WrywVVtVryc*5VyriV*y«?yr*V.-yryr**^ 

****    PURPOSE  -  ROUTINE  SUBTRACTS  MAT2  FROM  MAT1  AND  PUTS  THE 

*  RESULT  IN  A  THIRD  MATRIX. 

REAL  MAT1(10,10),MAT2(10,10),RMAT(10,10) 
INTEGER  IRR, IRC 

CALL  SMULTI(-1.0,MAT2,I2R,I2C) 

CALL  ADD(MAT1,I1R,I1C,MAT2,I2R,I2C,RMAT,IRR,IRC) 

IRR  =  IIR 

IRC  =  I1C 

RETURN 
END 

y?  yiyryryryryryryryryryry.-yrycyryTyryryryryTyTyryryjyryryTyTyTycyfyTyTyrycyryryfyryryryryryryfyryryTyTyrycyTyf 

SUBROUTINE   TPOSE(MATl , IIR, I1C,RMAT, IRR, IRC) 
Vr  yryryrycyryryryT>vyryryryryfyryTyryryryryTyfyryryryTyryryTyryryfyryrrvyryTy?yryryryfyryTyryTyryryTyTyr'5vyryr 

****   PURPOSE -SUBROUTINE  TRANSPOSES  A  MATRIX  AND  PUTS  THE  RESULT 
INTO  A  NEW  MATRIX 

REAL  MAT1(10,10),RMAT(10,10) 
INTEGER  I, J, IRR, IRC 


DO  93  1=1, I1C 

DO  94  J=1,I1R 

RMAT(I,J)  = 

MAT1(J,I) 

94 

CONTINUE 

93 

CONTINUE 
IRR  =  I1C 
IRC  =  IIR 
RETURN 

END 
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B.     MULTICHANNEL  ALGORITHM 


C 

c 
c 
c 
c 

* 

* 
y? 

c 

c 


yryryry?yry?ycyry?y?yfyry?yryryryrywcyryryryry?yTyr**yryTyryr* 

*  PAUL  DAL  SANTO   8/15/88 

*  TWO-CHANNEL  SYSTEM  IDENTIFICATION  ALGORITHM 

*  PROGRAM  CALCULATES  THE  AR  AND  MA  PARAMETERS 

*  BASED  UPON  THE  INPUT  AND  OUTPUT  DATA  OF  A 

*  TEST  SYSTEM. 

*  SHIFT  USED  TO  CALC  RYY  FOR  THE  LINEAR 


•y.- 
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* 
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* 
* 
* 

* 

* 

Vc 
* 

* 


* 

■if 


* 


*  OF  THE  AR  PARAMETERS  IS  READ  FROM  * 

*  THE  COEFF  DATA  FILE  * 

*  * 

*  NUMBER  OF  ITERATIONS  IS  READ  FROM  COEFF  * 

*  DATA  FILE  * 

*  * 

*  * 


ideMt 


RAVDAT 


VARIABLE  DEFINITIONS 


•kirie-ic 


ARRAY  WHICH  CONTAINS  THE  INPUT  AND  OUTPUT  DATA 

OF  THE  SYSTEM  UNDER  TEST 
El        ARRAY  FOR  STORING  THE  STOPPING  PARAMETER  AND  THE 

COEFFICIENT  ERROR 
ARRD      ARRAY  FOR  STORING  THE  AR  PARAMETER  ESTIMATES 
ARRN      ARRAY  FOR  STORING  THE  MA  PARAMETER  ESTIMATES 
RYYM      AUTOCORRELATION  MATRIX  OF  THE  OUTPUT  DATA 
RYUM,RUYM   CROSSCORRELATION  MATRICES 
RUUM      AUTOCORRELATION  MATRIX  OR  THE  INPUT  DATA 
RUUINV    INVERSE  OF  THE  AUTOCORRELATION  MATRIX  OF  THE  INPUT  DATA 
RYYINV    INVERSE  OF  THE  AUTOCORRELATION  MATRIX  OF  THE  OUTPT  DATA 
DM        VECTOR  OF  CURRENT  MA  PARAMETER  ESTIMATES 
CM        VECTOR  OF  CURRENT  AR  PARAMETER  ESTIMATES 
TRUNRM    FUNCTION  WHICH  CALCULATES  THE  NORM  OF  THE  TRUE  VALUES 

OF  THE  TEST  SYSTEM'S  PARAMETERS 
WKMAT     WORKING  MATRIX  FOR  DOING  MATRIX  INVERSIONS 
X  TPO     TRANSPOSE  OF  THE  VECTOR  OF  INPUT  DATA 
Y         MATRIX  FOR  THE  OUTPUT  OF  THE  TEST  SYSTEM 
U         INPUT  TO  THE  TEST  SYSTEM 

COEFM     VECTOR  OF  TRUE  COEFFICIENTS  OF  THE  TEST  SYSTEM 
ITERA     NUMBER  OF  ITERATIONS  TO  PERFORM 
CORRLN    LENGTH  OF  CORRELATIONS  TO  USE  TO  CALCULATE  RYY,  RUU 

RYU,  AND  RUY 
YSIZE     ORDER  OF  THE  AR  PART  OF  THE  TEST  SYSTEM 
USIZE     ORDER  OF  THE  MA  PART  OF  THE  TEST  SYSTEM 
KTIME     THE  CURRENT  ITERATION 
SHFT      SIZE  OF  THE  STARTING  LAG  FOR  LINEAR  PREDICTION  OF 

THE  AR  PARAMETERS 
SEED      INITIALIZATION  PARAMETER  FOR  IMSL  ROUTINE  WHICH 

GENERATES  RANDOM  GAUSSIAN  NUMBERS 

INTEGER  VARIABLES  THAT  END  WITH  R  CONTAIN  THE  ROW  SIZE  OF 

A  PARTICULAR  ARRAY.   INTEGER  VARIABLES  THAT  END  WITH  C  CONTAIN 

THE  COLUMN  SIZE  OF  A  PARTICULAR  ARRAY. 


»'..  - '»  «£*  mS* 


VARIABLE  DECLARATIONS 


ititMc 


COMMON  /D/  RAWDAT(2,0: 2000), El(2,0:  1000), 
+ARRD(0:  1000,5) ,ARRN(0:  1000,5) 

REAL  RYYM(5,5),RYUM(5,5),RUUM(5,5),RUYM(5,5),RUUINV(5,5) 
REAL  RYYINV(5,5),DM(5,5),CM(5,5),TRUNRM 

INTEGER  RYYR , RYYC , RUUR , RUUC , RYUR , RYUC , RUYR , RUYC 
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DUA0037 
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INTEGER  DR,DC,CR,CC 

REAL  WKMAT(12,12),X  TPO( 10 , 10) , Y( 2,2) ,U( 1) ,COEFM( 10 , 10) 
INTEGER  WKR , WKC , XTR , XTC , COEFR , COEFC , ERR 

INTEGER  ITERA,CORRLN,YSIZE,USIZE,KTIME,SHFT 

DOUBLE  PRECISION  SEED 

*  TEMPORATY  MATRICES  FOR  PERFORMING  CALCUALTIONS 

REAL  T1M(5,5),T2M(5,5),T3M(5,5) 
INTEGER  T1R,T1C,T2R,T2C,T3R,T3C 

*  BEGIN  MAIN  PROGRAM 

* 

*  READ  IN  THE  SIZE  OF  THE  TEST  SYSTEM,  THE  NUMBER  OF  ITERATIONS 

*  TO  PERFORM  AND  THE  BEGINNING  LAG  OF  LINEAR  PREDICTION  OF  THE 

*  AR  PARAMETERS 

READ(4,*,END=22)  YSIZE , USIZE , COEFR, COEFC , ITERA,SHFT 

*  READ  IN  THE  TRUE  VALUES  OF  THE  COEFFICIENTS  OF  THE  TEST  SYSTEM 

CALL  RDMAT(COEFM, COEFR, COEFC) 

*  INITIALIZE  ROW  AND  COLUMN  SIZES  OF  AS  WELL  AS  OTHER  VARIABLES 

CORRLN  =  500 

RUUR  =  USIZE 

RUUC  =  USIZE 

RYYR  =  YSIZE  +  1 

RYYC  =  YSIZE  +  1 

T3R  =  YSIZE 

T3C  =  YSIZE 

DR  =  USIZE 

DC  =  1 

CR  =  YSIZE  +  1 

CC  =  1 

XTR  =  COEFC 

XTC  =  COEFR 

SEED  =  888. 0D1 

II  =  0 

Y(l,l)  =  0.0 

E1(1,0)  =0.0 

DIVISR  =  TRUNRM(COEFM, COEFR, YSIZE, USIZE) 

*  ZERO  OUT  THE  CORRELATION  MATRICES,  THE  MA  PARAMETER  VECTOR  AND 

*  VECTOR  OF  INPUT  DATA 

CALL  INIT( RYYM , RYYR , RYYC ,0.0) 

CALL  INIT( RUUM, RUUR, RUUC, 0.  0) 

CALL  INIT(RYUM, RYYR, RUUC, 0.0) 

CALL  INIT(RUYM, RUUR, RYYC, 0.0) 
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CALL  INIT(DM,DR,DC,0.0) 
CALL  INIT(X  TPO,XTR,XTC,0.  0) 

*  INITIALIZE  THE  AR  PARAMETER  VECTOR 

CALL  INITD(CM,CR,CC,1.0) 
CALL  PRMAT(CM,CR,CC) 

*  RUN  THE  FILTER  FOR  2000  TIME  STEPS  TO  GENERATE  OUTPUT  DATA 
■h 

DO  70  KTIME  =  0,2000 


* 

* 
* 


70 
15 


GET  VALUE  FOR  U(K).   U  IS  A  GAUSSIAN  RANDOM  VARIABLE. 

CALL  GGNML  (SEED,1,U) 
SHIFT  U(K)  INTO  X  VECTOR  TO  CREATE  X(K) 

CALL  SHIFT(X  TPO,XTR,XTC,YSIZE,USIZE,Y( 1, 1) ,U( 1) ) 
CALCULATE  VALUE  OF  Y(K) 

CALL  MULTI(X  TPO,XTR,XTC ,COEFM,COEFR,COEFC , Y, 1 , 1) 

SAVE  THE  INPUT  AND  OUTPUT  DATA 

RAVDAT(1, KTIME)  =  Y(l,l) 
RAWDAT(2,KTIME)  =  U(l) 

CONTINUE 

FORMAT  (2X,I4,2X,F8.  5,2X,F8.5) 


*  CALCULATE  THE  CORRELATION  MATRICES 

CALL  AUTCOR( 1 ,RYYM, RYYR, RYYC, CORRLN) 
CALL  AUTCORC 2 , RUUM , RUUR , RUUC , CORRLN) 
CALL  CRSCOR( 1,2, RYUM , RYYR , RUUC , CORRLN ) 
CALL  CRSCOR( 2,1, RUYM , RUUR , RYYC , CORRLN) 

*  INVERT  THE  AUTOCORRELATION  MATRICES  OF  THE  OUTPUT  AND  INPUT  DATA 

CALL  LINV2F( RYYM , RYYR , RYYC , RYYINV , 0 , WKMAT , ERR ) 
CALL  LINV2F( RUUM , RUUR , RUUC , RUUINV , 0 , WKMAT , ERR ) 

*  MULTIPLY  THE  INVERSE  OF  THE  AUTOCORRELATION  MATRICES  BY  THEIR 

*  RESPECTIVE  CROSSCORRELATION  MATRICES 

CALL  MULTIC RUUINV, RUUR, RUUC, RUYM, RUUR, RYYC, TIM, T1R, TIC) 
CALL  MULTI(RYYINV, RYYR, RYYC, RYUM, RYYR, RUUC, T2M,T2R,T2C) 

*  ESTIMATE  THE  AR  PARAMETERS  BY  LINEAR  PREDICTION 

IF  (SHFT. GE.  1)  THEN 

CALL  C0RLA4(DR, CM, CR,CC,T3M, CORRLN, SHFT) 
ENDIF 

WRITE  (3,26)  II,DM(1,1),DM(2,1),DM(3,1),DM(4,1),DM(5,1) 
WRITE  (3,16)  II,CM(1,1),CM(2,1),CM(3,1),CM(4,1),CM(5,1) 
WRITE  (3,*)  ' 
CALL  SAVE(1,0,II,DM,DR,DC) 
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CALL  SAVE(2,0,II,CM,CR,CC) 

*  CALCULATE  THE  COEFFICIENT  ERROR 

CALL  ERR2(COEFM,COEFR,CM,CR,DM,DR,II,YSIZE,USIZE,DIVISR) 

*  BEGIN  THE  ITERATIVE  PROCEDURE 


DO  100  II  =  1,ITERA 

CALCULATE  THE  VALUES  OF  THE  PARAMETERS  AT  ITERATION  II 

CALL  MULTI(T1M,T1R,T1C,CM,CR,CC,DM,DR,DC) 
CALL  MULTI(T2M,T2R,T2C,DM,DR,DC,CM,CR,CC) 
CALL  SAVE( 1,0,11, DM, DR, DC) 
CALL  SAVE(2,0,II,CM,CR,CC) 

CALCULATE  THE  STOPPING  PARAMETER 
CALL  ERROR(CM,CR,CC,DM,DR,DC,II) 

CALCULATE  THE  COEFFICIENT  ERFOR 

CALL  ERR2(COEFM,COEFR,CM,CR,DM,DR,II,YSIZE,USIZE,DIVISR) 


26 
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WRITE  (3,26)  II,DM(1,1),DM(2,1),DM(3,1),DM(4,1),DM(5,1) 

FORMAT  (I4,1X,'NUM  COEF  =' ,5( IX, F9.  6) ) 

WRITE  (3,1')  II,CM(1,1),CM(2,1),CM(3,1),CM(4,1),CM(5,1) 

FORMAT  (I4,1X,'DNM  COEF  =' ,5( IX, F9.  6) ) 

WRITE  (3/0  '    ERROR  =  \E1(1,II),'   COEF  ERROR  =  ',E1(2,II) 

WRITE  (3,*)  ' 


CM(1,1)  =  1.  0 
100   CONTINUE 
*  END  OF  THE  ITERATIVE  PROCEDURE 


22 

* 
* 


PLOT  THE  PARAMETER  VALUES 

CALL  PL0T1(ITERA,DR,CR,C0EFM) 

STOP 
END 


91 


************************************************* 

FUNCTION  TRUNRM  (MAT1 ,M1R,YSZ ,USZ) 

*************************************  ************ 


REAL  MAT1(M1R,1) 
INTEGER  YSZ,USZ 

VALUE  =  0 

DO  91  I  =  1,YSZ  +  USZ 

VALUE  =  VALUE  +  (MAT1( I , 1) )**2 
CONTINUE 

TRUNRM  =  SQRT( VALUE  +  1) 
RETURN 
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END 

•it  ft'ieirMtiric'k'ic&irJeiciticific'kif&itiritic'it'fcititicieitie&irJticftJcir'kftiricJtic 

SUBROUTINE  AUTCOR( IFST,MAT1 ,M1R, MIC, CORRLN) 

*  THIS  SUBROUTINE  CALCULATES  THE  AUTOCORRELATION  MATRIX  USING 

*  CORRELATIONS  OF  SIZE  CORRLN 

COMMON  /D/  RAWDAT(2,0: 2000) ,E1(2,0: 1000), 
+ARRD(0: 1000,5) ,ARRN(0:  1000,5) 
REAL  MAT1(M1R,M1C) 
INTEGER  I, J, K, CORRLN 

CALL  INIT(MAT1,M1R,M1C,0.0) 


* 
* 

* 


CALC  CORRELATIONS  ALONG  THE  FIRST  ROW  OF  THE  MATRIX.   THE 

SHIFT  =  ABS( NUMBER  OF  THE  COLUMN  -  1) 

LENGTH  OF  CORR  =  D  LENGTH  +  ABS( 1  -  THE  NUMBER  OF  THE  COLUMN) 

ONCE  CORR  HAVE  BEEN  CALCULATED  FOR  THE  FIRST  ROW,  THEY  CAN 
BE  COPIED  INTO  OTHER  ROWS.   HORIZ  DISTANCE  OF  THE  PARTICULAR 
ELEMENT  FROM  THE  MAIN  DIAGONAL  DETERMINES  WHICH  CORR  TO 
COPY.  THIS  DISTANCE  IS  GIVEN  BY  ABS(ROW  NUMBER  -  COLUMN 
NUMBER). 

CORRELATIONS  START  200  POINTS  FROM  THE  BEGINNING  OF  THE 
DATA  TO  ELIMINATE  THE  TRANSIENT  OF  THE  TEST  SYSTEM. 


93 

90 


DO  90  J  =  1,M1C 

ENDPT  =  200  +  CORRLN  -  ABS( 1  -  J) 

DO  93  K  =  200, ENDPT 

MAT1(1,J)  =  MAT1(1,J)  +  RAWDAT(IFST,K-(1-J))*RAWDAT(IFST,K) 

CONTINUE 

MAT1(1,J)  =  MAT1(1,J)/(C0RRLN  +  1  -  ABS(1  -  J)) 
CONTINUE 


92 
91 


DO  91  I  =  2, MIR 

DO  92  J  =  1,M1C 

MAT1(I5J)  =  MAT1(1,1+ABS(I-J)) 

CONTINUE 
CONTINUE 


RETURN 
END 

SUBROUTINE  C0PY(MAT1 ,M1R,M1C,MAT2 ,M2R,M2C) 

*  THIS  SUBROUTINE  COPIES  A  MATRIX  INTO  A  SECOND  MATRIX. 

REAL  MAT1(M1R,M1C),MAT2(M2R,M2C) 
INTEGER  I, J 

DO  90  I  =  1,M1R 
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DO  900  J  =  1,M1C 

MAT2(I,J)  =  MAT1(I,J) 
900      CONTINUE 
90    CONTINUE 
M2R  =  MIR 
M2C  =  MIC 
RETURN 
END 

SUBROUTINE  CRSCOR( IFST, ISND,MAT1 ,M1R,M1C,C0RRLN) 

*  THIS  SUBROUTINE  CALCULATES  THE  CROSS CORRELATION  MATRIX  USING 

*  CORRELATIONS  OF  LENGTH  CORRLN 

COMMON  /D/  RAWDAT(2,0:  2000), El(2,0:  1000), 
+ARRD(0: 1000,5) ,ARRN(0:  1000,5) 

REAL  MAT1(M1R,M1C) 
INTEGER  I, J, K, CORRLN 

CALL  INIT(MAT1,M1R,M1C,0.  0) 

*  FOR  CROSS  CORRELATION  MUST  CALC  EACH  ELEMENT  OF  THE  CORR  MATRIX 
SEPARATELY.   CORRELATIONS  START  200  POINTS  FROM  THE  BEGINNING 

*  OF  THE  DATA  TO  ELIMINATE  THE  TRANSIENT  OF  THE  TEST  SYSTEM. 

DO  94  I  =  1,M1R 
DO  95  J  =  1,M1C 

ENDPT  =  200  +  CORRLN  -  ABS(I  -  J) 
IF  (J.GE.  I)  THEN 

DO  96  K  =  200, ENDPT 

MAT1(I,J)  =  MAT1(I,J)  +  RAWDAT(IFST,K-(I-J)) 
+*RAWDAT(ISND,K) 

96  CONTINUE 
ELSE 

DO  97  K  =  200, ENDPT 

MAT1(I,J)  =  MAT1(I,J)  +  RAWDAT(IFST,K)* 
+RAWDAT( I SND , K+( I -J) ) 

97  CONTINUE 
ENDIF 

MAT1(I,J)  =  MAT1(I,J)/(C0RRLN  +  1  -  ABS(I-J)) 
95       CONTINUE 
94    CONTINUE 

RETURN 

END 

SUBROUTINE  C0RLA4( ISIZE,MAT2 ,M2R,M2C,MAT3, CORRLN, SHIFTT) 

*  THIS  SUBROUTINE  CALCULATES  THE  CORRELATION  MATRIX 

*  AND  THE  CORRELATION  VECTOR  USED  FOR  LINEAR  PREDICTION  OF  THE 

*  AR  PARAMETERS.   IT  THEN  CALCULATES  THE  INITIAL  ESTIMATE  OF  THE 

*  AR  PARAMETERS  AND  PASSES  THEM  BACK  TO  THE  MAIN  PROGRAM  IN  MAT2. 
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COMMON  /D/  RAWDAT( 2,0: 2000), El( 2,0: 1000), 
+ARRD(0: 1000 ,5) , ARRN( 0: 1000,5) 
REAL  MAT2(M2R,M2C),MAT3(M2R-1,M2R-1) 
REAL  T2M( 5,1), T3M( 5,1), T4M( 5,5), WKMAT( 9,9) 
INTEGER  ERR,T1R,T1C,T2R,T2C,T3R,T3C,ISIZE,SHIFTT,C0RRLN 

*  GENERATE  THE  FIRST  ROW  OF  THE  RYY  MATRIX. 

*  SHIFTS  ARE  GREATER  THAN  THE  ORDER  OF  THE 

*  NUMERATOR. 

M3R  =  M2R-1 

M3C  =  M3R 

CALL  INIT( MAT3 , M3R , M3C , 0. 0 ) 

DO  90  J  =  1,M3C 

ENDPT  =  200  +  CORRLN  -  SHIFTT 
DO  91  K  =  200, ENDPT 

MAT3(1,J)  =  MAT3(1,J)  +  RAWDAT( 1 ,K+SHIFTT)*RAWDAT( 1 ,K) 

91  CONTINUE 

MAT3(1,J)  =  MAT3(l,J)/(CORRLN-SHIFTT+l) 
SHIFTT  =  SHIFTT  +  1 
90    CONTINUE 

*  COPY  ELEMENTS  FROM  THE  FIRST  ROW  INTO  OTHER  LOCATIONS 

DO  92  I  =  2,M3R 
DO  93  J  =  1,M3C 

MAT3(I,J)  =  MAT3(1,1+ABS(I-J)) 

93  CONTINUE 

92  CONTINUE 

*  GENERATE  THE  RYY  VECTOR  BY  COPYING  ELEMENTS  OF  THE  RYY  MATRIX 

T2R  =  M3C 

T2C  =  1 

CALL  FILL( 1 ,T2R- 1,1, 1 ,T2M ,T2R ,1,2, 1 ,MAT3 ,M3R ,M3C) 

*  GENERATE  THE  LAST  ELEMENT  IN  THE  RYY  CORRELATION  VECTOR. 

FINELE  =  0.0 

ENDPT  =  200  +  CORRLN  -  SHIFTT 

DO  94  K  =  200, ENDPT 

FINELE  =  FINELE  +  RAWDAT( 1,K+SHIFTT)*RAWDAT( 1,K) 

94  CONTINUE 

FINELE  =  FINELE/ (CORRLN+1 -SHIFTT) 

*  COPY  THE  FINAL  ELEMENT  INTO  THE  VECTOR 
T2M(T2R,1)  =  FINELE 

*  CALCULATE  THE  INITIAL  ESTIMATE  OF  THE  AR  PARAMETERS 

CALL  LINV2F(MAT3,M3R,M3C,T4M,0,WKMAT,ERR) 

CALL  MULTI ( T4M , M3R , M3C , T2M , T2R , T2C , T3M , T3R , T3C ) 

*  COPY  THE  AR  PARAMETERS  INTO  THE  RETURN  ARGUMENT 
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* 
* 

95 


* 
* 


MAT2(1,1)   =   1.  0 
DO   95   L  =   1,3 

MAT2(L  +   1,1)   =  T3M(L,1) 

MAT2(L,1)   =  T3M(L-1,1) 

WRITE   (3,*)   MAT2(L,1) 
CONTINUE 

RETURN 
END 

SUBROUTINE  ERROR(MATl ,M1R,M1C,MAT2,M2R,M2C , ITNUM) 


*  THIS   SUBROUTINE   CALCULATES   THE   STOPPING   PARAMETER. 

COMMON   /D/   RAWDAT(2,0: 2000), El( 2,0:  1000), 
+ARRD(0:  1000 ,5) ,ARRN(0:  1000,5) 
REAL  MAT1(M1R,M1C) ,MAT2(M2R,M2C) 
INTEGER   I,J,K 

ERVAL  =0.0 
XVAL  =  0.  0 
ZVAL  =0.0 


DO  90   I   =  400,450 

DO  91   J  =  MIR, 1,-1 

XVAL  =  XVAL  +  MAT1(J,1)*RAWDAT(1,I+M1R-J) 
CONTINUE 
DO  92   K  =  M2R,1,-1 

ZVAL  =  ZVAL  +  MAT2(K,1)*RAWDAT(2,I+M2R-K) 
CONTINUE 


91 


92 


90 


ERVAL  =  ERVAL  +  (XVAL-ZVAL)**2 
XVAL  =  0.  0 
ZVAL  =   0.0 
CONTINUE 

E 1(1, ITNUM)   =  ERVAL/ 51 


RETURN 
END 

»•.  ju  jl  j-  j.  Jm  jl  .'-  ju  j*  ju  -' -  -J-  ■»<-  j-  -'-  -^  -/-  •.*-  -*-  -v  J-  y-  -j-  j-  «  y-  -J-  -'-  y-  y-  J-  -•-  y-  y-  yr  y-  -*r  yr  y r  y-  J  -  -y  -/"':  v  r 

SUBROUTINE  ERR2( MAT 1,M1R,MAT2,M2R,MAT3,M3R, ITNUM, 
+YSZ,USZ,DIVISR) 

*  THIS   SUBROUTINE  CALCUALTES  THE  COEFFICIENT  ERROR. 

COMMON  /D/   RAWDAT( 2,0:  2000), El( 2,0:  1000), 
+ARRD(0: 1000,5) ,ARRN(0:  1000,5) 
REAL  MAT1(M1R,1),MAT2(M2R,1),MAT3(M3R,1) 
INTEGER   I,YSZ,USZ 

ERRVAL  =  (1    -   MAT2(1,1))**2 
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90 


DO  90  I  =  1,YSZ 

ERRVAL  =  ERRVAL  +  (MAT1(I,1)  +  MAT2( 1+1 , 1) )**2 
CONTINUE 


DO  92  I  =  1,USZ 

ERRVAL  =  ERRVAL  +  (MAT1( I+YSZ, 1)  -  MAT3(I,1))**2 
92    CONTINUE 

E1(2,ITNUM)  =  SQRT( ERRVAL) /DIVISR 

RETURN 
END 

SUBROUTINE  FILL( I , J,K,L,MAT1 ,M1R,M1C,R2 ,C2,MAT2 ,M2R,M2C) 

Vr       V!rifcyfVlryc:fc:V*:feyT*A*<IWr^yJryr*y*Vr^ 

*  THIS  ROUTINE  FILLS  MAT1  FROM  MAT2.   POSITIONS 

*  IN  MAT1  FROM  (I,K)  TO  (J,L)  ARE  FILLED  WITH  AN  EQUAL  NUMBER 

*  OF  ELEMENTS  FROM  MAT2  STARTING  AT  POSITION  (R2,C2). 

REAL  MAT1(M1R,M1C) ,MAT2(M2R,M2C) 
INTEGER  ROW,COL,R2,C2,C22 

C22  =  C2 


91 


90 


DO  90  ROW  =  I, J 
DO  91  COL  =  K,L 

MATl(ROW,COL)  =  MAT2(R2,C2) 
C2  =  C2  +  1 
CONTINUE 
R2  =  R2  +  1 
C2  =  C22 
CONTINUE 


RETURN 
END 

SUBROUTINE  LIM2( EMAX , ESTP , CEMAX , CESTP , ITERA) 

■JU  .  ■ ,  JL  Jg  J+  J •  . '..  JL  J-  JU  JL  *'v-  JU  JU  J>  JL  J-  JU  JL  J-  JL  .'.  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  JL  J-  JL  JL  JL  JL  JL  J-  J-  J.  JL  JL  JL  V-  JL  J-  JL  J-  J- 

*  ROUTINE  CALCULATES  THE  LIMITS  FOR  THE  GRAPH  OF  THE  STOPPING 

*  PARAMETER  AND  THE  COEFFICIENT  ERROR. 

*  STOPPING  PARAMETER  LIMITS  ARE  RETURNED  IN  EMAX  AND  ESTP. 

*  COEFFICIENT  ERROR  LIMITS  ARE  RETURNED  IN  CEMAX  AND  CESTP. 

COMMON  /D/  RAWDAT(2,0:  2000), El(2,0:  1000), 
+ARRD(0:  1000,5) ,ARRN(0:  1000,5) 
REAL  EMAX, ESTP 
INTEGER  ITERA 

EMAX  =  0.  0 

DO  90  I  =  0, ITERA 

IF  (El( 1,1). GT. EMAX)  THEN 
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EMAX  =  E1(1,I) 
END  IF 

90  CONTINUE 

EMAX  =  1.  25  *  EMAX 
ESTP  =  EMAX/5 
CEMAX  =  0.  0 

DO  91  I  =  0, ITERA 

IF  (E1(2,I).GT.  CEMAX)  THEN 

CEMAX  =  E1(2,I) 
ENDIF 

91  CONTINUE 

CEMAX  =  1. 25  *  CEMAX 
CESTP  =  CEMAX/5 
RETURN 
END 

SUBROUTINE  PL0T1( ITERA, ICURVN, ICURVD,MAT1) 

Vc        Vr?VVriVVr?VV-VrVfVfVrV-V-V-'5C'5ViV-V'jV'!V-)'--V'5V'!VVr'5VVc-.'c-.V-,V 

*  THIS  ROUTINE  GENERATES  SEPARATE  PLOTS  OF  THE  MA 

*  AND  AR  PARAMETERS.   IT  THEN  REPLOTS  THE  THESE  CURVES  ALONG 

*  WITH  THE  STOPPING  PARAMETER.   FINALLY  IT  PLOTS  THE  STOPPING 

*  PARAMETER  AND  THE  COEFFICIENT  ERROR  ON  THE  SAME  GRAPH. 

COMMON  /D/  RAWDAT(2,0: 2000), El(2,0:  1000), 
+ARRD(0: 1000 ,5) ,ARRN(0: 1000,5) 
REAL  X(0: 1000), Y(0:  1000) ,MAT1( 10 , 10) 
REAL  STP,RNMIN,RNSTEP,RNMAX,RITERA 
INTEGER  I , J , I 1R , I 1C , VAL , ITERA , ICURV 

CALL  LIMITS( ICURVN, ICURVD,DMAX,DMIN,DSTEP, 
+NMAX,NMIN,NSTEP, ITERA) 
CALL  LIM2( EMAX, ESTP, CEMAX, CESTP, ITERA) 

*  GENERATE  THE  ITERATION  NUMBER 


DO  90  I  =  0, ITERA 
X(I)  =  I 
Y(I)  =0.0 

CONTINUE 


90 


*  CALCULATE  X  AXIS  LABELING  INTERVAL 

STP  =  ITERA/ 10. 

*  SET  THE  DEVICE  CALLS  FOR  ALL  OF  THE  PLOTS 

CALL  SHERPA('MCGRAF   ','A',3) 
CALL  RESET('ALL') 

*  SECTION  1 

*  THIS  SECTION  GRAPHS  THE  NUMERATOR  AND  DENOMINATOR 

*  PARAMETERS  ON  SEPARATE  GRAPHS  ON  THE  SAME  PAGE. 


DUA07190 
DUA07200 
DUA07220 

DUA07230 
DUA07240 
DUA07280 

DUA07290 
DUA07310 
DUA07320 
DUA07330 
DUA07350 

DUA07360 
DUA07370 
DUA07380 
DUA07390 
DUA07400 
DUA07670 
DUA07680 
DUA07690 
DUA07700 
DUA07710 
DUA07720 
DUA07730 
DUA07740 
DUA07750 
DUA07760 
DUA07770 
DUA07780 
DUA07790 
DUA07800 
DUA07810 
DUA07820 
DUA07830 
DUA07840 
DUA07870 
DUA07880 
DUA07890 
DUA07900 
DUA07910 
DUA07920 
DUA07930 
DUA07940 
DUA07950 
DUA07970 
DUA07980 
DUA07990 
DUA08000 
DUA08030 
DUA08050 
DUA08060 
DUA08070 
DUA08080 
DUA08090 
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*  SET  UP  THE  PLOT  FOR  THE  NUMERATOR  COEFF 

CALL  PAGE(  8.  5,11.0) 

CALL  HEIGHT(0.  14) 

CALL  HVROT('AUTO') 

CALL  XINTAX 

CALL  PHYSOR(  1.5,6.  0) 

CALL  AREA2D(5.0,3.5) 

CALL  COMPLX 

CALL  YAXANG(O) 

CALL  XNAME( ' ITERATIONS? ' , 100) 

CALL  YNAMEC  PARAMETER  VALUES?  '  ,  100) 

CALL  MESSAG(' (A)?' ,100,2.4,-0. 8) 

CALL  THKFRM(0. 03) 

CALL  FRAME 

CALL  GRAF(0. ,STP, ITERA,NMIN,NSTEP,NMAX) 

*  PLOT  THE  NUMERATOR  PARAMETERS 

DO  93  J  =  1,ICURVN 
DO  94  I  =  0,ITERA 

Y(I)  =  ARRN(I,J) 
94       CONTINUE 

CALL  CURVE (X,Y, ITERA,0) 
93    CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUE  OF  THE  PARAMETERS 

CALL  DASH 

DO  9  7  K  =  ICURVD,ICURVD+ICURVN-1 
DO  98  J  =  0,ITERA 
Y(J)  =  MAT1(K,1) 
98       CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
97    CONTINUE 

CALL  ENDGR(O) 

*  SET  UP  THE  PLOT  FOR  THE  DENOMINATOR  PARAMETERS 

CALL  RESET('DASH') 

CALL  PHYSORC  1.5,1.  5) 

CALL  AREA2D(5.  0,3.5) 

CALL  XNAME( ' ITERATIONS? ' , 100) 

CALL  YNAMEC 'PARAMETER  VALUES? ', 100) 

CALL  MESSAG('(B)?' ,100,2.4,-0. 8) 

CALL  THKFRM(0.03) 

CALL  FRAME 

CALL  GRAF(0. ,STP, ITERA,DMIN,DSTEP,DMAX) 

*  PLOT  THE  DENOMINATOR  PARAMETERS 

DO  95  J  =  1,ICURVD 
DO  96  I  =  0,ITERA 

Y(I)  =  ARRD(I,J) 
96       CONTINUE 


DUA0810C 

DUA0811C 

DUA0812C 

DUA0813C 

DUA0814C 

DUA0815C 

DUA0816C 

DUA0817G 

DUA0820C 

DUA0821C 

DUA0824C 

DUA0825C 

DUA0826C 

DUA0827C 

DUA0828C 

DUA0829C 

DUA0832C 

DUA0834C 

DUA0835C 

DUA0836C 

DUA0837C 

DUA0838C 

DUA0839C 

DUA0840C 

DUA0841C 

DUA0842C 

DUA0843C 

DUA0844C 

DUA0845C 

DUA0846C 

DUA0847C 

DUA08480 

DUA08490 

DUA0850Q 

DUA08510 

DUA08520 

DUA08530 

DUA08540 

DUA08550 

DUA08560 

DUA08570 

DUA08590 

DUA08600 

DUA08620 

DUA08630 

DUA08640 

DUA08650 

DUA08660 

DUA08680 

DUA08700 

DUA08710 

DUA08720 

DUA08730 

DUA08740 

DUA08750 

DUA08760 
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95 


CALL  CURVE (X,Y, ITERA,0) 
CONTINUE 


*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUES  OF  THE  DENOM  PARAMETERS 


CALL  DASH 

DO  99  K  =  1,ICURVD-1 

DO  100  J  =  0,ITERA 
Y(J)  =  -MAT1(K,1) 

CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
CONTINUE 


100 
99 

105 


DO  105  J  =  0,ITERA 
Y(J)  =1.0 

CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 

CALL  ENDPL(O) 


*  SECTION  2 

*  THIS  SECTION  PUTS  THE  STOPPING  PARAMETER  ON  THE 

*  GRAPHS  OF  THE  NUMERATOR  AND  DENOMINATOR  PARAMETERS. 

*  SET  UP  THE  PLOT  FOR  THE  NUMERATOR  PARAMETERS 

CALL  RESET('DASH') 

CALL  HWROT('AUTO') 

CALL  XINTAX 

CALL  PHYSOR(  1.5,6.0) 

CALL  AREA2D(5.  0,3.5) 

CALL  COMPLX 

CALL  YAXANG(O) 

CALL  XNAME( ' ITERATIONS$ ' , 100) 

CALL  YNAME(' PARAMETER  VALUES$ ' , 100) 

CALL  MESSAG( '(A)$' ,100,2.  4,-0.  8) 

CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF(0. ,STP,ITERA,NMIN,NSTEP,NMAX) 

*  PLOT  THE  NUMERATOR  PARAMETERS 

DO  200  J  =  1,ICURVN 
DO  201  I  =  0,ITERA 

Y(I)  =  ARRN(I,J) 
201       CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
200    CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUES  OF  THE  PARAMETERS 

CALL  DASH 

DO  202  K  =  ICURVD,ICURVD+ICURVN-1 
DO  203  J  =  O.ITERA 
Y(J)  =  MAT1(K,1) 

203       CONTINUE 


DUA08770 
DUA08780 
DUA08790 
DUA08800 
DUA08810 
DUA08820 
DUA08830 
DUA08840 
DUA08850 
DUA08860 
DUA08870 
DUA08880 
DUA08890 
DUA08900 
DUA08910 
DUA08920 
DUA08930 
DUA08940 
DUA08950 
DUA08960 
DUA08970 
DUA08980 
DUA08990 
DUA09010 
DUA09020 
DUA09030 
DUA09040 
DUA09050 
DUA09060 
DUA09090 
DUA09100 
DUA09130 
DUA09140 
DUA09150 
DUA09160 
DUA09170 
DUA09180 
DUA09190 
DUA09210 
DUA09220 
DUA09230 
DUA09240 
DUA09250 
DUA09260 
DUA09270 
DUA09280 
DUA09290 
DUA09300 
DUA09310 
DUA09320 
DUA09330 
DUA09340 
DUA09350 
DUA09360 
DUA09370 
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202 


CALL  CURVE(X,Y,ITERA,0) 
CONTINUE 


*  PLOT  THE  STOPPING  PARAMETER  ON  THE  SAME  GRAPH 

CALL  DOT 

CALL  YGRAXS(0.0,ESTP,EMAX, 3.  5 ,' STOPPING  PARAMETER? ' , 
+  -100,5.0,0.  0) 

DO  204  J  =  0,ITERA 
Y(J)  =  E1(1,J) 
204   CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
CALL  ENDGR(O) 


*  SET  UP  THE  PLOT  FOR  THE  DENOMINATOR  PARAMETERS 

CALL  RESET('DOT') 

CALL  PHYS0R(1.  5,1.5) 

CALL  AREA2D(5.  0,3.5) 

CALL  XNAME( ' ITERATIONS $ ' , 100) 

CALL  YNAME(' PARAMETER  VALUES? ', 100) 

CALL  MESSAG( '(B)?' ,100,2.4,-0.8) 

CALL  THKFRM(0.03) 

CALL  FRAME 

CALL  GRAF(0. ,STP, ITERA,DMIN,DSTEP,DMAX) 

*  PLOT  THE  DENOMINATOR  PARAMETERS 

DO  205  J  =  1,ICURVD 
DO  206  I  =  0,ITERA 

Y(I)  =  ARRD(I,J) 
206      CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
205   CONTINUE 

*  PLOT  DASHED  LINES  FOR  THE  TRUE  VALUES  OF  THE  PARAMETERS 


CALL  DASH 

DO  207  K  =  1,ICURVD-1 
DO  208  J  =  0,ITERA 

Y(J)  =  -MAT1(K,1) 
CONTINUE 

CALL  CURVE (X,Y,ITERA,0) 
CONTINUE 


208 
207 

209 


DO  209  J  =  0,ITERA 

Y(J)  =  1.0 
CONTINUE 
CALL  CURVE(X,Y,ITERA,0) 


*  PLOT  THE  STOPPING  PARAMETER  ON  THE  SAME  GRAPH 
CALL  DOT 


DUA09380 
DUA09390 
DUA09400 
DUA09410 
DUA09420 
DUA09430 
DUA09440 
DUA09450 
DUA09460 
DUA09470 
DUA09480 
DUA09490 
DUA09500 
DUA09510 
DUA09520 
DUA09530 
DUA09540 
DUA09550 
DUA09560 
DUA09570 
DUA09580 
DUA09590 
DUA09610 
DUA09620 
DUA09630 
DUA09640 
DUA09650 
DUA09660 
DUA09680 
DUA09690 
DUA09700 
DUA09  710 
DUA09720 
DUA09730 
DUA09740 
DUA09750 
DUA09760 
DUA09770 
DUA09780 
DUA09790 
DUA09800 
DUA09810 
DUA09820 
DUA09830 
DUA09840 
DUA09850 
DUA09860 
DUA09870 
DUA09880 
DUA09890 
DUA09900 
DUA09910 
DUA09920 
DUA09930 
DUA09940 
DUA09950 
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210 


CALL  YGRAXS(0.  0 ,ESTP,EMAX, 3. 5 ,* STOPPING  PARAMETER$', 
+-100,5. 0,0. 0) 

DO  210  J  =  0,ITERA 

Y(J)  =  El(l, J) 
CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 
CALL  ENDPL(O) 


*  SECTION  3 

*  THIS  SECTION  PLOTS  THE  STOPPING  PARAMETER  AND  THE  COEFFICIENT 

*  ERROR  ON  THE  SAME  GRAPH. 

*  SETUP  THE  PLOT  FOR  THE  STOPPING  PARAMETER 

CALL  DOT 

CALL  HWROT('AUTO') 

CALL  PHYSOR(  1.5,6.  0) 

CALL  AREA2D(5.0,3.5) 

CALL  XNAME( ' ITERATIONS$ ' , 100) 

CALL  YNAME( ' STOPPING  PARAMETER? ' , 100) 

CALL  THKFRM(0.  03) 

CALL  FRAME 

CALL  GRAF(0.  , STP, ITERA,0.  ,ESTP,EMAX) 


DO  306  J  =  0,ITERA 
Y(J)  =  E1(1,J) 
306   CONTINUE 

CALL  CURVE (X,Y,ITERA,0) 


*  PLOT  THE  COEFFICIENT  ERROR  ON  THE  SAME  GRAPH 

CALL  CHNDOT 

CALL  YGRAXS(0.0,CESTP,CEMAX, 3.  5,' COEFFICIENT  ERROR?', 
+-100,5. 0,0. 0) 

DO  307  J  =  0,ITERA 
Y(J)  =  E1(2,J) 
307   CONTINUE 

CALL  CURVE(X,Y,ITERA,0) 

CALL  ENDPL(O) 

CALL  DONEPL 

RETURN 

END 

SUBROUTINE  SAVE(VAL,M,K,MAT1,M1R,M1C) 

*  ROUTINE  SAVES  PARAMETER  ESTIMATES  IN  EITHER  ARRD  OR  ARRN 

*  DEPENDING  UPON  THE  VALUE  OF  VAL. 

COMMON  /D/  RAWDAT(2,0:  2000), El(2,0:  1000), 
+ARRD(0: 1000,5) ,ARRN(0:  1000,5) 


DUA09960 
DUA099  70 
DUA09980 
DUA09990 
DUA10000 
DUAinoiO 
DUA10020 
DUA 10030 
DUA10040 
DUA10050 
DUA10060 
DUA10070 
DUA10080 
DUA10090 
DUA10100 
DUA10110 
DUA10120 
DUA10130 
DUA10140 
DUA10150 
DUA10160 
DUA10180 
DUA10190 
DUA10200 
DUA10220 
DUA10230 
DUA10240 
DUA10250 
DUA10260 
DUA10270 
DUA10280 
DUA10290 
DUA10300 
DUA10310 
DUA10320 
DUA10330 
DUA10340 
DUA 10 350 
DUA10360 
DUA10370 
DUA10380 
DUA10390 
DUA10400 
DUA10410 
DUA10420 
DUA10430 
DUA10440 
DUA10860 
DUA10870 
DUA10880 
DUA10890 
DUA10900 
DUA10910 
DUA10920 
DUA10930 
DUA10940 
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REAL  MAT1(M1R,M1C) 
INTEGER   I,J,K,M,VAL 


DO   90    I   =   1,M1R 

DO   900   J  =    1,M1C 

IF   (VAL.EQ.  1)   THEN 

ARRN(K.I)   =  MAT1(I,J) 

ELSEIF   (VAL.EQ.  2)   THEN 

ARRD(K,I)   =  MAT1(I,J) 

ENDIF 

900 

CONTINUE 

90 

CONTINUE 

RETURN 

END 

DUA1096C 
DUA1097C 
DUA1098C 
DUA1099C 
DUA1100C 
DUA1101C 
DUA1102C 
DUA1103C 
DUA1104C 
DUA1105C 
DUA1106C 
DUA1107C 
DUA1108C 
DUA1109C 
DUA1110C 


C.     SUBPROGRAMS  COMMON  TO  BOTH  SYSTEM  IDENTIFICATION 
ALGORITHMS  SUB00010 


SUB00020 


* 


920 
92 


* 


95 

94 


SUBROUTINE   ADD   (MAT1 , IR1 , IC1 ,MAT2 , IR2 , IC2 ,RMAT, IRR, IRC) 

jl  y  -  jl  y,  -  ■  -  jl  y.  jl  jl  -J- J-  j>  y-  jl  y  -  yf  jl  j*  jl  jl .  ■ .  .'-  jl  -i-  jl  jl  j-  J-  jl  jl  JL jl  jl  jl  y-  y  -  j  -  jl  JL  jl  y -  J-  jl  y-  y  -  y  -  y-  y  - jl 


SUB0004C 
SUB0005C 
SUB0006C 
SUB0007C 
THIS   SUBROUTINE  ADDS  TWO  EQUAL  SIZE  MATRICIES  AND  PUTS  THE  RESULT  SUB0008C 


IN  A  THIRD  MATRIX. 

REAL  MAT1(IR1,IC1),MAT2(IR2,IC2),RMAT(IR1,IC1) 
INTEGER   I, J, IRR, IRC 

DO  92   1=1, IR1 

DO   920  J=1,IC1 

RMAT(I,J)   =  MAT1(I,J)   +  MAT2(I,J) 

CONTINUE 
CONTINUE 
IRR  =   IR1 
IRC   =   IC1 

RETURN 
END 


- ■ .  JL JL  .." -  JL  JL  JL  - ' .  JL.  JL.  JL  JL  JL 


SUBROUTINE   1NIT(MAT1 , MIR, MIC , INITVL) 


THIS   SUBOUTINE   INITIALIZES  A  MATRIX  TO   INITVL 


REAL  MATKM1R, MIC),  INITVL 
INTEGER   I, J 

DO  94   1=1, MIR 

DO  95   J=1,M1C 

MAT1(I,J)=INITVL 

CONTINUE 
CONTINUE 
RETURN 


SUB0009C 

SUB0010C 

SUB0013C 

SUB0014C 

SUB0016C 

SUB0017C 

SUB0018C 

SUB0019C 

SUB0020C 

SUB0021C 

SUB0022C 

SUB0023C 

SUB0024C 

SUB00250 

SUB00260 

SUB00280 

SUB00290 

SUB00300. 

SUB00310 

SUB00320 

SUB00330 

SUB00340 

SUB00350 

SUB00360 

SUB00380 

SUB00390 

SUB00400 

SUB00410 

SUB00420 

SUB00430 
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* 

* 


95 

94 

* 
* 

* 


91 

90 


END 

SUBROUTINE  INITD(MAT1 , MIR, MIC, INITVL) 

THIS  SUBROUTINE  INITIALIZES  A  MATRIX  TO  AS  A  DIAGONAL  MATRIX 
WHOSE  DIAGONAL  ELEMENTS  EQUAL  INITVL.  ' 

REAL  MAT1(M1R, MIC), INITVL 
INTEGER  I, J 

DO  94  1=1, MIR 
DO  95  J=1,M1C 

IF  (I.EQ.  J)  THEN 

MAT1(I,J)=INITVL 
ELSE 

MAT1(I,J)=0.0 
ENDIF 
CONTINUE 
CONTINUE 
RETURN 
END 

SUBROUTINE  LIMITS(NSZ,DSZ,DMAX,DMIN,DSTEP,NMAX,NMIN,NSTEP, ITERA) 

ROUTINE  CALCULATES  THE  LIMITS  FOR  THE  GRAPHS 
CALCULATES  DENOMINATOR  AND  NUMERATOR  LIMITS  SEPARATELY 
IN  PREPARATION  FOR  MAKING  TWO  GRAPHS 

COMMON  /D/  RAWDAT(2,0: 2000), El(2,0: 1000), 
+ARRD(0:  1000 ,5) , ARRN( 0:  1000,5) 
REAL  DMAX,DMIN,DSTEP,NMAX,NMIN,NSTEP 
INTEGER  DSZ,NSZ 

CALCULATE  THE  DENOMINATOR  LIMITS 

DMAX  =1.0 
DMIN  =  0.  0 

DO  90  I  =  1,DSZ 

DO  91  J  =  0, ITERA 

IF  ((ARRD(J,I)).GT.  DMAX)  THEN 
DMAX  =  ARRD(J,I) 

ENDIF 

IF  ((ARRD(J,I)).LT.  DMIN)  THEN 
DMIN  =  ARRD(J,I) 

ENDIF 

CONTINUE 
CONTINUE 

IF  (DMAX. GT.  0)  THEN 

DMAX  =  1.25  *  DMAX 
ELSE 


SUB00440 
SUB00450 
SUB00470 
SUB00480 
SUB00490 
SUB00500 
SUB00510 
SUB00520 
SUB00530 
SUB00540 
SUB00550 
SUB00570 
SUB00580 
SUB00590 
SUB00600 
SUB00610 
SUB00620 
SUB00630 
SUB00640 
SUB00650 
SUB00660 
SUB00670 
SUB00680 
SUB00690 
SUB00710 
SUB00720 
SUB00730 
SUB00740 
SUB00750 
SUB00760 
SUB00770 
SUB00780 
SUB00790 
SUB00800 
SUB00820 
SUB00830 
SUB00840 
SUB00850 
SUB00860 
SUB00870 
SUB00880 
SUB00890 
SUB00900 
SUB00910 
SUB00920 
SUB00930 
SUB00940 
SUB00950 
SUB00960 
SUB00970 
SUB00980 
SUB00990 
SUB01000 
SUB01010 
SUB01020 
SUB01030 
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* 


DMAX  =  0.  0 
ENDIF 

IF  (DMIN. GT.  0)  THEN 

DMIN  =  0.0 
ELSE 

DMIN  =  1. 25  *  DMIN 
ENDIF 
DSTEP  =  (DMAX  -  DMIN)/5 

CALCULATE  THE  NUMERATOR  LIMITS 

NMAX  =  0.  0 
NMIN  =  0.  0 
DO  92  I  =  1,NSZ 

DO  93  J  =  0,ITERA 

IF  (ARRN(J,I).GT.  NMAX)  THEN 
NMAX  =  ARRN(J,I) 

ENDIF 


IF  (ARRN(J,I).LT.  NMIN)  THEN 

NMIN  =  ARRN(J,I) 

ENDIF 

93 

CONTINUE 

92 

CONTINUE 

IF  (NMAX. GT. 0)  THEN 

NMAX  =  1.  25  *  NMAX 

ELSE 

NMAX  =  0.  0 

ENDIF 

IF  (NMIN.  GT.  0)  THEN 

NMIN  =  0. 0 
ELSE 

NMIN  =  1.  25  *  NMIN 
ENDIF 
NSTEP  =  ABS(NMAX  -  NMIN)/5 

RETURN 
END 

SUBROUTINE  MULTI  (MAT1 ,M1R,M1C,MAT2,M2R}M2C,RMAT,M3R,M3C) 

ROUTINE  MULTIPLIES  TWO  MATRICES  AND  PUT  THE  RESULT  IN  A 
THIRD  MATRIX. 

REAL  MAT1(M1R,M1C) ,MAT2(M2R,M2C) ,RMAT(M1R,M2C) 
INTEGER  I^KjIRRJRC 

CALL  INIT(RMAT,M1R,M2C,0.  0) 

DO  91  1=1, MIR 

DO  910  J=1,M2C 


SUB01040 

SUB01050 

SUB0106C 

SUB01070 

SUB01080 

SUB01090 

SUB01100 

SUB0111C 

SUB01130 

SUB01140 

SUB01150 

SUB01160 

SUB0117C 

SUB0118C 

SUB01190 

SUB01200 

SUB01210 

SUB01220 

SUB01230 

SUB01240 

SUB01250 

SUB01260 

SUB01270 

SUB01280 

SUB01290 

SUB01300 

SUB01310 

SUB01320 

SUB01330 

SUB01340 

SUB01350 

SUB01360 

SUB01370 

SUB01380 

SUB01390 

SUB01400 

SUB0141O 

SUB0143O 

SUB01440 

SUB0145O 

SUB0146O 

SUB01470: 

SUB01490 

SUB01500 

SUB01510: 

SUB0152O 

SUB01530 

SUB01540 

SUB01550i 

SUB015801 

SUB01590 

SUB01600 

SUB01610 

SUB01620 

SUB01630 

SUB01640 
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DO  9100  K=1,M1C 

RMAT( I , J)=RMAT( I , J)  +MAT1( I ,K)*MAT2(K, J) 
9100        CONTINUE 
910      CONTINUE 
91    CONTINUE 
M3R  =  MIR 
M3C  =  M2C 
RETURN 
END 

"k  *iV*}V*}ViV}VV«W«V*iVVr^V}V?V:fc*iV**yfiV?V}VitoV 

SUBROUTINE  PRMAT(MAT1 , I1R, I1C) 


SUB01650 
SUB01660 
SUB01670 
SUB01680 
SUB01690 
SUB01700 
SUB01710 
SUB01720 
SUB01730 
SUB01740 
IVA04910 
IVA04920 
IVA04930 


*  SUBROUTINE  PRINTS  A  MATRIX  OUT  TO  THE  FILE  DEFINED 

*  AS  UNIT  3 


REAL  MAT1(10,10) 
INTEGER  I, J 

DO  92  I  =  1,I1R 

WRITE(3,302)  (MAT1(I,J),J  =  1,I1C) 
302         FORMAT  (7(2X,F8.5)) 
92       CONTINUE 
RETURN 
END 

SUBROUTINE  RDMAT(MAT1 , MIR, MIC) 

*  ROUTINE  READS  A  MATRIX  FROM  FILE  SPECIFIED  AS  UNIT  4. 

REAL  MAT1(M1R,M1C) 
INTEGER  I, J 

*  READ  IN  MATRIX 

DO  92  I  =  1,M1R 

READ  (4,*)  (MAT1(I,J),J=1,M1C) 
C301         FORMAT(10F3.  1) 
92       CONTINUE 
RETURN 
END 

SUBROUTINE  SHIFT(MAT1 ,M1R,M1C3DSIZE,NSIZE ,0UTDAT, INDAT) 

*  ROUTINE  SHIFTS  NEW  INPUT  AND  OUTPUT  VALUES  INTO  DATA  VECTOR 

*  OF  THE  TEST  SYSTEM.   THE  OLDEST  VALUES  ARE  LOST  TO  MAKE 

*  ROOM  FOR  THE  NEW  VALUES. 

REAL  MAT1(M1R,M1C) ,OUTDAT( 1) , INDAT( 1) 
INTEGER  J,DSIZE,NSIZE,NSTART 


IVA04970 
1VA04980 
IVA04990 
IVA05010 
IVA05020 
IVA05030 
IVA05040 
IVA05050 
IVA05060 
IVA05070 
SUB01760 
SUB01770 
SUB01780 
SUB01790 
SUB01800 
SUB01810 
SUB01840 
SUB01850 
SUB01860 
SUB01870 
SUB01880 
SUB01890 
SUB01900 
SUB01910 
SUB01920 
SUB01930 
SUB01940 
SUB01950 
SUB01970 
SUB01980 
SUB01990 
SUB02000 
SUB02010 
SUB02020 
SUB02030 
SUB02060 
SUB02070 
SUB02080 
SUB02090 
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NSTART  =  DSIZE  +  1 
DO  92  J  =  DSIZE, 2,-1 

MAT1(1,J)  =  MAT1(1,J-1) 

92  CONTINUE 

MAT1( 1,1)  =  OUTDAT(l) 

DO  93  J  =  M1C,NSTART+1,-1 
MAT1(1,J)  =  MAT1(1,J-1) 

93  CONTINUE 

MAT 1(1, NSTART)  =  INDAT(l) 
WRITE(12,19)  (MAT1(1,J),J=1,M1C) 
19    FORMAT  (7(1X,F8.4)) 

RETURN 
END 

it  ititititititititititititititititititititititititititititititirititititititititititit 

SUBROUTINE  SMULTI  ( CONST, MAT1 , MIR, MIC) 

•ft  it  sV  it  jV  «V  "V  ;Y  sV  itit  it  ■>  c  A  -.V  ic  ft  >V  »V  A  A  iV  A  ~V  -V  it  it  it  it  it  it  it  it  it  it  it  itit  it  it  it  it  it 

*     ROUTINE  MULTIPLIES  A  MATRIX  BY  A  SCALAR. 

REAL  MAT1(M1R,M1C), CONST 
INTEGER  I, J 

DO  93  1=1, MIR 

DO  930  J=1,M1C 

MAT1(I,J)  =  MAT1(I,J)  *  CONST 
930      CONTINUE 
93    CONTINUE 

RETURN 

END 

SUB02480 


SUB02100 

SUB02110 

SUB02120 

SUB02130 

SUB02140 

SUB02150 

SUB02160 

SUB02170I 

SUB02180 

SUB02190 

SUB02200 

SUB02210 

SUB02220 

SUB02230 

SUB02240 

SUB02260 

SUB02270 

SUB02280 

SUB02290 

SUB02300 

SUB02310 

SUB02320 

SUB02350 

SUB02370 

SUB02380 

SUB02390 

SUB02400 

SUB02410 

SUB02420 

SUB02430 

SUB02440 

SUB02450 

SUB02460 
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